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Abstract 

The density-matrix renormalization group method (DMRG) has established itself over the last 
decade as the leading method for the simulation of the statics and dynamics of one-dimensional 
strongly correlated quantum lattice systems. In the further development of the method, the re- 
alization that DMRG operates on a highly interesting class of quantum states, so-called matrix 
product states (MPS), has allowed a much deeper understanding of the inner structure of the 
DMRG method, its further potential and its limitations. In this paper, I want to give a detailed 
exposition of current DMRG thinking in the MPS language in order to make the advisable im- 
plementation of the family of DMRG algorithms in exclusively MPS terms transparent. I then 
move on to discuss some directions of potentially fruitful further algorithmic development: while 
DMRG is a very mature method by now, I still see potential for further improvements, as exem- 
plified by a number of recently introduced algorithms. 
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1. Introduction 

Strongly correlated quantum systems on low-dimensional lattices continue to pose some of 
the most interesting challenges of modern quantum many-body physics. In condensed matter 
physics, correlation effects dominate in quantum spin chains and ladders, in frustrated magnets 
in one and two spatial dimensions, and in high-temperature superconductors, to name but a few 
physical systems of interest. More recently, the advent of highly controlled and tunable strongly 
interacting ultracold atom gases in optical lattices has added an entirely new direction to this 
fieldfl']. 

Both analytically and numerically, these systems are hard to study: only in very few cases 
exact analytical solutions, for example by the Bethe ansatz in one dimension, are available|0,|3t]. 
Perturbation theory fails in the presence of strong interactions. Other approaches, such as field 
theoretical approaches, have given deep insights, for example regarding the Haldane gap physics 
of integer-spin antiferromagnetic chains 1 4], but make potentially severe approximations that 
must ultimately be controlled by numerical methods. Such algorithms include exact diagonal- 
ization, quantum Monte Carlo, series expansions or coupled cluster methods. 

Since its invention in 1992 by Steve Whitejsl j6t], the density-matrix renormalization group 
(DMRG) has firmly established itself as the currently most powerful numerical method in the 
study of one-dimensional quantum lattices|0, After initial studies of the static properties 
(energy, order parameters, «-point correlation functions) of low-lying eigenstates, in particu- 
lar ground states, of strongly correlated Hamiltonians such as the Heisenberg, t-J and Hubbard 
models, the method was extended to the study of dynamic properties of eigenstates, such as dy- 
namical structure functions or frequencv-dependent conductivitiesllgj IToi \ \, T2l. At the same 
time, its extension to the analysis of two-dimensional classical |T3| and one-dimensional quan- 
tum 1IT4I [15I1 transfer matrices has given access to highly precise finite-temperature information 
on classical two-dimensional and quantum one-dimensional systems; more recently the transfer 
matrix variant of DMRG has also been extended to dynamics at finite temperature 1 16). It has 
even been extended to the numerically much more demanding study of non-Hermitian (pseudo-) 
Hamiltonians emerging in the analysis of the relaxation towards classical steady states in one- 
dimensional systems far from equiUbriumlB H [H, Q . 
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In many applications of DMRG, the accuracy of results is essentially limited only by machine 
precision, even for modest numerical resources used, quite independent of the detailed nature of 
the Hamiltonian. It is therefore not surprising that, beyond the extension of the algorithm to 
more and more problem classes, people wondered about the physical origins of the excellent 
performance of DMRG and also whether the success story could be extended to the study of 
real-time dynamics or of two-dimensional systems. 

In fact, both questions are intimately related: as was realized quite soon, DMRG is only mod- 
erately successful when applied to two-dimensional lattices: while relatively small systems can 
be studied with high accuracv ll2ll l22i I23I l24i l25i I26I, l27i I28I1 . the amount of numerical resources 
needed essentially increases exponentially with system size, making large lattices inaccessible. 
The totally different behaviour of DMRG in one and two dimensions is, as it turned out, closely 
related 1291 [30I1 to the different scaling of quantum entanglement in many-body states in one and 
two dimensions, dictated by the so-called area laws (for a recent review, see 13111 ). 

In this paper, I will stay within the framework of one-dimensional physics; while the gen- 
erahzations of DMRG to higher dimensions reduce naturally to DMRG in one dimension, the 
emerging structures are so much richer than in one dimension that they are beyond the scope of 
this work. 

In an originally unrelated development, so-called matrix product states (MPS) were discov- 
ered as an interesting class of quantum states for analytical studies. In fact, the structure is so 
simple but powerful that it is no surprise that they have been introduced and used under a variety 
of names over the last fifty or more years (most notably perhaps by Baxter i32il ). In the present 
context, the most relevant prehistory is arguably given by the exact expression of the seminal 
one-dimensional AKLT state in this form [33,,34,.,35.1 . which gave rise to extensive studies of the 
translationally invariant subclass of MPS known as finitely correlated states lf36ll . This form was 
then subsequently used in a variety of contexts for analytical variational calculations, e.g. for 
spin-1 Heisenberg antiferromagnets[37..38, 39, 40] and ferrimagnets 14 U 14211 . 

The connection between MPS and DMRG was made in two steps. In a first step, Ostlund 
and Rommer realized that the block-growth step of the so-called infinite-system DMRG 
could be expressed as a matrix in the form it takes in an MPS. As in homogeneous systems this 
block-growth step leads to a fixed point in the thermodynamic limit, they took the fixed point 
matrix as building block for a translationally invariant MPS. In a further step, it was recognized 
that the more important finite-system DMRG leads to quantum states in MPS form, over which it 
variationally optimizes ll4411 . It was also recognized that in traditional DMRG the state class over 
which is variationally optimized changes as the algorithm progresses, such that if one demands 
in some sense "perfect" variational character, a small change to the algorithm is needed, which 
however was found to increase (solvable) metastability problems ll45[ l46ll . 

It remains a curious historical fact that only a few of the DMRG practicioners took this 
development very seriously up to about 2004 when Cirac, Verstraete, Vidal and coworkers started 
to explore the power of MPS very systematically. While it was considered useful for conceptual 
purposes, surprisingly little thought was given to rethinking and reexpressing real-life DMRG 
implementations purely in the MPS language; arguably, because the overwhelming majority of 
conventional DMRG applications (i.e. ground states for quantum chains with open boundary 
conditions) hardly profits. What was overlooked is that it easily opens up the way to powerful 
extensions to DMRG hard to see and express in conventional DMRG lan g ua g e. 

A non-exhaustive list of extensions would list real-time evolutions 
54], also at finite temperature ['sT, , the efficient use of periodic boundary conditionsr56, 57|, 
5811 . reliable single-site DMRG Ii461 . numerical renormalization group (NRG) applications 16211 . 
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infinite-system algorithmsll63. 64, 65], continuous systemsf66'l, not talking at all about progress 
made in higher dimensions starting with 16711 using a generalization of the MPS state class L68,l . 

The goal of this paper cannot be to provide a full review of DMRG since 1992 as seen from 
the perspective of 2010, in particular given the review|7], which tries to provide a fairly extensive 
account of DMRG as of summer 2004. I rather want to limit myself to more recent developments 
and phrase them entirely in the language of matrix product states, focussing rather on the nuts 
and bolts of the methods than showing a lot of applications. My hope would be that this review 
would allow newcomers to the field to be able to produce their own codes quickly and get a 
firm grasp of the key building blocks of MPS algorithms. It has overlaps with the introductions 
1 6^ 7^ in the key methods presented, but focuses on different extensions, some of which arose 
after these pap ers, and in many places tries to be more explicit. It takes a different point of 
view than 117 III , the first comprehensive exposition of DMRG in 1998, because at that time the 
connection to MPS (though known) and in particular to quantum information was still in many 
ways unexploited, which is the focus here. Nevertheless, in a first "historical" step, I want to 
remind readers of the original way of introducing DMRG, which does not make use of the idea 
of matrix product states. This should make older literature easily accessible, but one can jump to 
Section|4]right away, if one is not interested in that. 

In a second step, I wiU show that any quantum state can be written exactly in a very specific 
form which is given by the matrix product states akeady alluded to. In fact, the restriction to 
one dimension will come from the fact that only in this case MPS are numerically manageable. 
I will highlight special canonical forms of MPS and establish their connection to the singular 
value decomposition (SVD) as a mathematical tool and the Schmidt decomposition as a compact 
representation of quantum states. After this I will explain how MPS are a natural framework 
for decimation schemes in one dimension as they occur in schemes such as DMRG and Wil- 
son's NRG. As a simple, but non-trivial example, I will discuss the AKLT state in its MPS form 
explicitly. We then move on to discuss explicitly operations with MPS: overlaps, normaliza- 
tion, operator matrix elements, expectation values and MPS addition. These are operations one 
would do with any quantum state; more MPS-specific are methods for bringing them into the 
computationally convenient canonical forms and for approximating an MPS by another one of 
smaller dimension. I conclude this exposition of MPS with discussing the relationship and the 
conversions between the MPS notation I favour here, an alternative notation due to Vidal, and the 
DMRG way of writing states; this relatively technical section should serve to make the literature 
more accessible to the reader 

The MPS ideas generalize from states to the representation of operators, so I move on to 
discuss the use of matrix product operators (MPO)[51, 70, 72, TS^ 74]. As far as I can see, all 
operators of interest to us (ranging from local operators through bond evolution operators to full 
Hamiltonians) find a very compact and transparent formulation in terms of MPO. This leads to 
a much cleaner and sometimes even numerically more accurate formulation of DMRG-related 
algorithms, but their usage is not yet very widely spread. 

Admittedly, at this point the patience of the reader may have been stretched quite a bit, as 
no real-world algorithm e.g. for ground state searches or time evolutions has been formulated 
in MPS language yet; but it will become obvious that a lot of cumbersome numerical details of 
DMRG algorithms have been hidden away neatly in the MPS and MPO structures. 

I will discuss ground state algorithms, discussing the equivalences and differences between 
DMRG with one or two center sites and fully MPS-based algorithms, including improvements 
to avoid trapping. I will focus on finite systems with open boundary conditions, where these 
methods excel. 
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After this, I move on to time-dependent methods for dynamics, for pure and mixed states. 
After a discussion of the basic algorithms and their subtle differences, I will focus on the key 
problem of extending the time-range of such simulations: The possibility to calculate highly 
accurate real-time and imaginary-time evolutions of complex quantum many-body states has 
been particularly exciting for many people, also because it arrived just in time for studying highly 
tunable ultracold atom systems. While this development has already led to numerous interesting 
insights and applications, it was quickly realized that the time-range of time-dependent DMRG 
and related methods is limited by entanglement growth in quantum states out of equilibrium, 
such that long-time physics is out of reach. In this context, interesting progress in trying to go 
beyond has been achieved recently. 

The review concludes with two further axes of development. I will start out by discussing 
the connection between DMRG and Wilson's NRG, showing how NRG can be expressed in a 
very concise fashion as well as be improved in various directions. This closes an interesting 
historical loop, as the utter failure of NRG for homogeneous one-dimensional quantum lattices 
as opposed to quantum impurity models mapped to special non-homogeneous one-dimensional 
quantum lattices was at the starting point of White's invention of DMRGITS*]. 

I continue by looking at infinite-size algorithms using MPS that work directly in the thermo- 
dynamic limit, one based on time evolution (iTEBD)||63|]. The other (iDMRG)165] is an exten- 
sion of infinite-system DMRG algorithm, which has had an interesting history: in many early 
discussions of DMRG it was presented as the key aspect of DMRG, with finite-system DMRG 
as a practitioners' add-on to further improve numerical accuracy. Later, it was recognized that 
applying finite-system DMRG is essential even for qualitative correctness in many cases, and 
infinite-system DMRG was seen as just a warm-up procedure. Only recently, McCulloch [65ll 
pointed out a way how to turn infinite-system DMRG into a highly efficient tool for producing 
thermodynamic limit states for homogeneous systems. 

Last but not least, I will give an outlook on further applications of MPS that I could not cover 
here. 

2. Density -matrix renormalization group (DMRG) 

2.1. Infinite-system and finite-system algorithms 

As a toy model, let us consider an (anisotropic) S - ^ Heisenberg antiferromagnetic (/ = 1) 
spin chain of length L in one spatial dimension with external magnetic field h. 



We consider open boundary conditions (Fig.[Tl), which is well worth emphasizing: analytically, 
periodic boundary conditions are usually more convenient; many numerical methods do not re- 
ally depend strongly on boundary conditions, and some, like exact diagonalization, even become 
more eflicient for periodic boundary conditions. DMRG, on the other hand, prefers open bound- 
ary conditions. 

It should also be emphasized that, DMRG being a variational method in a certain state class, 
it does not suffer from anything like the fermionic sign problem, and can be applied to bosonic 
and fermionic systems alike. 

The starting point of DMRG is to ask for the ground state and ground state energy of H. We 
can ask this question for the thermodynamic limit L — > oo or more modestly for finite L. In the 




L-l 



L 



(1) 
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Figure 1 : Our toy model: a chain of length L with open ends, where a spin^ 
nearest neighbours. 



sits on each site and interacts with its 



first case, the answer is provided by infinite-system DMRG albeit with quite hmited precision; in 
the second case, an answer can be read off from infinite-system DMRG, but it is more adequate to 
run a two-step procedure, starting with infinite-system DMRG and continuing with finite-system 
DMRG. 

In any case, the numerical stumbling block is provided by the exponential growth of the 
Hilbert space dimension, in our example as d^, where d - lis the local state space dimension of 
a spin-i. 

2.2. Infinite-system DMRG 

Infinite-system DMRG deals with this problem by considering a chain of increasing length, 
usually L = 2, 4, 6, . . ., and discarding a sufficient number of states to keep Hilbert space size 
manageable. This decimation procedure is key to the success of the algorithm: we assume that 
there exists a reduced state space which can describe the relevant physics and that we can develop 
a procedure to identify it. The first assumption is typical for all variational methods, and we will 
see that indeed we are lucky in one dimension: for all short-ranged Hamiltonians in ID there is 
such a reduced state space that contains the relevant physics ! 

How is it found? In infinite-system DMRG (Fig.|2|, the buildup is carried out as follows: we 
introduce left and right blocks A and B, which in a first step may consist of one spin (or site) 
each, such that total chain length is 2. Longer chains are now built iteratively from the left and 
right end, by inserting pairs of spins between the blocks, such that the chain grows to length 4, 
6, and so on; at each step, previous spins are absorbed into the left and right blocks, such that 
block sizes grow as 1,2, 3, and so on, leading to exponential growth of the dimension of the 
fuU block state space as 2^, where { is the current block size. Our chains then always have a 
block-site-site-block structure, A»»B. 

Let us assume that our numerical resources are sufficient to deal with a reduced block state 
space of dimension D, where in practice D will be of (9(100) to (9(1000), and that for a block A 
of length { we have an effective description of block A in a D-dimensional reduced Hilbert space 
with orthonormal basis For very small blocks, the basis dimension may be less than D, 

for larger blocks some truncation must have occurred, which we take for granted at the moment. 
Let me mention right now that in the literature, D - which will turn out to be a key number - 
comes under a variety of names: in traditional DMRG literature, it is usually referred to as m (or 
M); more recent matrix product state literature knows both D and;^f. 

Within this framework, we may now ask, (i) what is the ground state of the current chain of 
length 2{ -H 2, also referred to as superblock, and (ii) how can we find the reduced Hilbert space 
of dimension D for the new blocks A» and "B. 
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Any state of the superblock A»»B can be described by 



\^)^ X '/'fl,lcr,icrs<jJfl>A|0")A|0")B|fl)B = ^lAMjJOAly^B, (2) 

aAO-AtTBaB iaJb 

where the states of the site next to A are in the set {|cr)^) of local state space dimension d, 
and analogously those of the site next to B. By numerical diagonalization we find the li/r) that 
minimizes the energy 

^ ^ {l/f\HA..B\<('} 

with respect to the Hamiltonian of the superblock, answering (i). To this purpose, we need some 
iterative sparse matrix eigensolver such as provided by the Lanczos or Jacobi-Davidson methods. 
Given that our Hamiltonian [Eq. ([TJ] is available in the full tensor product space, this minimiza- 
tion of course assumes that it can be expressed readily in the superblock basis. Let us postpone 
this question, as it is intimately related to growing blocks and decimating their state space. As 
the matrix dimension is d^D~, for typical D the eigenproblem is too large for direct solution, 
even assuming the use of quantum symmetries. As most matrix elements of short-ranged Hamil- 
tonians are zero, the matrix is sparse, such that the basic matrix-vector multiplication of iterative 
sparse matrix eigensolvers can be implemented very efficiently. We will discuss this question 
also further below. 

If we now take states {|/)^) as the basis of the next larger left block A», the basis dimension 
grows to dD. To avoid exponential growth, we truncate this basis back to D states, using the 
following procedure for block A» and similarly for "B: we consider the reduced density operator 
for A», namely 

Pa. = Tr.B|(/r)<(A| (pa.)»' = ^ ^ij<P*'j- (4) 

j 

The eigensystem of pA. is determined by exact diagonalization; the choice is to retain as reduced 
basis those D orthonormal eigenstates that have the largest associated eigenvalues w. If we call 
them \b)A, the vector entries are simply the expansion coefficients in the previous block-site basis, 
A{a(T\b)A- 

After an (approximate) transformation of all desired operators on A» into the new basis, the 
system size can be increased again, until the final desired size is reached. B is grown at the same 
time, for reflection-symmetric systems by simple mirroring. 

The motivation of the truncation procedure is twofold. The first one, which is on a weaker 
basis, is that we are interested in the states of A» contributing most to the ground state for A» em- 
bedded in the final, much larger, even infinitely large system. We approximate this final system 
ground state, which we don't know yet, to the best of our knowledge by that of A»»B, the largest 
superblock we can efliciently form. In that sense, we are bootstrapping, and the finite-system 
DMRG will, as we will see, take care of the large approximations this possibly incurs. Let me 
remark right now that in the MPS formulation of infinite-system DMRG a much clearer picture 
will emerge. The second motivation for the choice of states is on much firmer foundation: the 
above prescription follows both from statistical physics arguments or from demanding that the 
2-norm distance |||^^) - |iA)ti-uncll2 between the current ground state li/*) and its projection onto trun- 
cated block bases of dimension D, |iA)tmnc' is minimal. The most transparent proof follows from a 
singular value decomposition of the matrix with matrix elements ^ (aAa-A),{a-BaB) formed from the 
wave function coefficients ^'a^a-.^a-Bai,- The overall success of DMRG rests on the observation that 
even for moderate D (often only a few 100) the total weight of the truncated eigenvalues, given 



by the truncation error e = 1 - Yja>D w'a if we assume descending order of the Wg, is extremely 
close to 0, say 10"'" or less. 

Algorithmically, we may of course also fix a (small) e we are willing to accept and at each 
truncation choose D sufficiently large as to meet this requirement. Computational resources are 
used more efficiently (typically, we will have a somewhat smaller D towards the ends of chains 
because of reduced quantum fluctuations), the programming effijrt to obtain this flexibility is of 
course somewhat higher. 

An important aspect of improving the performance of any quantum algorithm is the exploita- 
tion of quantum symmetries, ranging from discrete symmetries like a Z2 mirror symmetry if 
ffie Hamiltonian is invariant under mirroring from left to right through Abelian symmetries like 
ffie t/(l) symmetry of particle number ("charge") or magnetization conservation to non- Abelian 
symmetries like ffie SU{2) symmetry of rotational invariance frMf^fl^ . A huge range of 



ffiese symmetries have been used successfully in numerous applications, with the most common 
ones being the two f/Y 1 ^symmetries of charge and magnetization conservation, but also SU{2) 
symmetries! 8^ 81 , 82, 83 1. Let us focus on magnetization and assume ffiat ffie total magnetiza- 



tion, M - commutes with the Hamiltonian, [H,M] - 0, such that eigenstates of H can 

be chosen to be eigenstates of M. Let us assume in particular that the good quantum number 
of the ground state is M = 0. If the block and site states are eigenstates of magnetization, then 
i^aAO-Aa-BOB + only if -^(|fl)A) + -^(|o")a) + A/(|o")b) + b) = 0; M is a short-hand for the 
magnetization of the respective blocks and sites, assuming that the states have magnetization as a 
good quantum number. This constraint allows to exclude a large number of coefficients from the 
calculation, leading to a huge speedup of the calculation and (less important) savings in memory. 

The decisive point is that if the states and |cr)^ are eigenstates of magnetization, so will 
be the eigenstates of the reduced density operator pA., which in turn will be the states |fl)^ of 
ffie enlarged block. As the local site states can be chosen to be such eigenstates and as the first 
block in the growth process consists of one local site, the eigenstate property then propagates 
through the entire algorithm. To prove ffie claim, we consider (pa.),,' - YjJ The states 

li}^ and are eigenstates of magnetization by construction, hence M(\i}^) + M{\j)g) = = 
M{\i')A) + M{\j)g) or M(|/)^) = M(|/')^). The density matrix therefore decomposes into blocks 
that are formed from states of equal magnetization and can be diagonalized block by block within 
these sets, such that its eigenstates are also eigenstates of magnetization. 

In practical implementations, the use of such good quantum numbers will be done by arrang- 
ing matrix representations of operators into block structures labelled by good quantum numbers, 
which are then combined to satisfy local and global constraints. While this is conceptually easy, 
ffie coding becomes more complex and will not be laid out here explicitly; hints will be given in 
the MPS sections. 

So far, we have postponed the question of expressing operators acting on blocks in the current 
block bases, in order to construct the Hamiltonian and observables of interest. Let us consider 
an operator O acting on site I, wiffi matrix elements (9°"' °"f = {cr(\Oi\a-'f). Without loss of gener- 
ality, assume that site £ is on the left-hand side and added into block A. Its construction is then 
initialized when block A grows from ^ - 1 — > ^, as 

{a(\a(-x(T ()(sT c\0\(T'^{at-ia-'y^. (5) 
Here, \ac) and are the effective basis states of blocks of length t and I - \ respectively. Of 
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course, updates are necessary during the further growth steps of block A, e.g. 



(flf+i|(9|fl^^j) = ^ <a^+i|fl^o-c+i)<a^|0|a^)<a'^cr^+i|flJ,^j). (6) 
It is important to realize that the sum in Eq. (|6]) must be split as 



^ <fl^+i|fl;o-f+i> ^<flf|(9|flP<fl;Cr(.+i|a^^j) 



(7) 



reducing the calculational load from O(D'^d) to 20(D^d). 

In Hamiltonians, operator products OP occur. It is tempting, but due to the many truncation 
steps highly imprecise, to insert the quantum mechanical one in the current block basis and to 
write 

{at\dP\a'f} = 2 {ae\6\ae}{a(\P\a'f}. NO! (8) 
iif 

The correct way is to update the first operator (counting from the left for a left block A) until the 
site of the second operator is reached, and to incorporate it as 

{a{\dP\a'i,) = ^ {a{\a(-icr[){ac-i\d\a'^_^){o-(\P\o-'i,){a'f_^o-'i,\a'f). (9) 

fl,_i,a;_i,iTf,a-'f 

The further updates are then as for single operators. Obviously, we are looking at a straightfor- 
ward sequence of (reduced) basis transformations, with a somewhat cumbersome notation, but it 
will be interesting to see how much simpler these formulae will look in the MPS language, albeit 
identical in content. 

The ultimate evaluation of expectation values is given at the end of the growth procedure as 
{(f/\d\4r) = ^ {il/\aAa-AO-BaB){aA\6\a'^}{a'^o-AcrBaB\if'}, (10) 

aAa'J^a■J^a■BaB 

where suitable bracketing turns this into an operation of order O(D^d^). An important special 
case is given if we are looking for the expectation value of a local operator that acts on one of 
the free sites • in the final block-site configuration. Then 

{i//\d\i//} ^ ^ {tff\aAO-AO-BaB}{o-A\d\o-'^}{aAcr'^o-BaB\i//}, (11) 

aACTACT'^^a-BaB 

which is an expression of order O(D^d^) (for many operators, even O(D^d^)). Given that D » d 
in all practical applications, such evaluations are computationally highly advantageous. 

Let me conclude these more technical remarks by observing that for an efficient construction 
of H from such operators, it is essential never to build it as a full matrix, but to make use of 
the specific block-site-site-block structure. Assume, for example, a term which contains one 
operator acting on A and one on B (this is in fact the most complicated case), h = OaOb- Then 



{aAO-AcrBaB\h\>/^} = ^ <flA|OA|a^) 



{aB\dB\a'g){a'^o-AO-Ba'g\ilf) 



(12) 
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Figure 2: The left and right half of the figure present the iterations taken in the infinite-system and finite-system DMRG 
procedures respectively. In both cases, new blocks are formed from integrating a site into a block, with a state space 
truncation according to the density-matrix prescription of DMRG. Whereas in the infinite- system version this growth 
happens on both sides of the chain, leading to chain growth, in the finite-system algorithm it happens only for one side 
at the expense of the other, leading to constant chain length. 



which is a sequence of two O(D^cfi) multiplications (instead of one naive O(D'^cfi) calculation) 
for all coefficients (fl^cr^cr^a^l/ilt^) and similarly for all other terms. 

If we stop infinite-system DMRG at some superblock size L, we can interpret the final wave- 
function in two ways. We can take it as an approximation to the exact state for the superblock 
of size L and evaluate expectation values. The accuracy is limited not only by the truncations, 
but also by the fact that the first truncations were carried out for extremely small superblocks: 
the choice of relevant short block states is likely to be a not too good approximation to those one 
would have chosen for these short blocks embedded in the final system of length L. 

Alternatively, we may ignore these boundary effects and focus on the central sites as an 
approximation to the behaviour of an infinitely large system, provided L is large enough and 
careful extrapolations of results to L — > oo are done. This is not so simple to do in a very con- 
trolled fashion, but as we will see, infinite-system DMRG can be used to build highly controlled 
translationally invariant (modulo, say, a unit cell of length 2) thermodynamic limit states. 



2.3. Finite-system DMRG 

Once the desired final system size is reached by infinite-system DMRG, it is important in 
all but the most trivial applications to follow up on it by the so-called finite-system DMRG 
procedure. This will not merely lead to some slight quantitative improvements of our results, 
but may change them completely: consider 18411 for an example where even the central physical 
statement changes: for a t-J-V-V model on a ladder with moderate hole-doping S = 0.1, an 
infinite-system DMRG calculation indicates the existence of alternately circulating currents on 
plaquettes that are triggered by an infinitesimal current at the left end of the ladder, a signal 
of a so-called t/-density wave state. Only after applying the finite-system algorithm it becomes 
obvious that this current is in fact exponentially decaying into the bulk, excluding this type of 
order (Fig.[3]|. 

The finite-system algorithm corrects the choices made for reduced bases in the context of a 
superblock that was not the system of interest (of final length L), but some sort of smaller proxy 
for it. 

What the finite-system algorithm does is the following (Figure |2|: it continues the growth 
process of (say) block B following the same prescription as before: finding the ground state of the 
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Figure 3: Plaquette current on a t-J-V-V -ladder, J = OAt, V = 3t, V = t (V nearest, V next-nearest neighbour 
interaction) at liole doping 6 = 0.1, system size 2 X 60, as induced by a finite boundary current on rung 1 . The absolute 
current strength is shown; whereas infinite-system DMRG and the first sweep indicate the generation of a long-ranged 
pattern, a fully converged calculation (here after 6 to 7 sweeps) reveals an exponential decay into the bulk. Taken from 
Ref. dl. 

superblock system, determining the reduced density operator, finding the eigensystem, retaining 
the D highest weight eigenstates for the next larger block. But it does so at the expense of block 
A, which shrinks (i.e. old shorter blocks A are reused). This is continued until A is so small as to 
have a complete Hilbert space, i.e. of dimension not exceeding D (one may also continue until A 
is merely one site long; results are not affected). Then the growth direction is reversed: A grows 
at the expense of B, including new ground state determinations and basis choices for A, until B 
is small enough to have a complete Hilbert space, which leads to yet another reversal of growth 
direction. 

This sweeping through the system is continued until energy (or, more precisely, the wave 
function) converges. The intuitive motivation for this (in practice highly successful) procedure 
is that after each sweep, blocks A or B are determined in the presence of an ever improved 
embedding. 

In practice, this algorithm involves a lot of book-keeping, as all the operators we need have to 
be maintained in the current effective bases which will change from step to step. This means that 
the truncated basis transformations determined have to be carried out after each step; operator 
representations in all bases have to be stored, as they will be needed again for the shrinking block. 

Another important feature is that for finding the ground state for each A»»B configuration 
one employs some iterative large sparse matrix eigensolver based on sequential applications of 
H on some initial starting vector To speed up this most time-consuming part of the algorithm, 
it is highly desirable to have a good prediction for a starting vector, i.e. as close as possible to 
the ultimate solution. This can be achieved by (approximately) transforming the result of the 
last step into the shifted A»»B configuration 12111 by applying two basis transformations: e.g. 
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A» and for a sweep to the right. The explicit formulae (see iQIzH]) can be derived 

by writing 

llA) = ^ ^aiaMac,2b,Jae)A\cr{+\)\a-{+2)\b{+2)B, (13) 

where |flf)^ and \bc+2)B ^6 the block states for block A comprising sites 1 through £ and block 
B comprising sites i + 3 through L (the label of the block states is taken from the label of 
the bond their ends cut; see Fig. [T3]) . and inserting twice an approximate identity, namely / - 

2flf+, \at+\)AA{ac+\\ and/ = Tjauibr^j. |o-m>l^f+3>B B<^f+3Ko"f+3l- One then obtains 

\^) = ^ lAflf+icrf+,trf,3Vj"f+l>Al0"^+2>|0"f+3>l^f+3>ij, (14) 

with 

{ae+\\afcre+\){bMcrM\bi+2)- (15) 

a((J"f+l6(+2 

The basis transformations required in the last equation are all available from previous steps in the 
DMRG procedure. A similar operation can be carried out for a sweep to the left. As we will see, 
this astute step, which led to drastic improvements in DMRG performance, is already implicit if 
one rewrites DMRG in the MPS language, such that we will not discuss it here at length. 

An important observation is that both the infinite-system and finite-system algorithm can also 
be carried out by inserting only a single explicit site •, hence one would study superblocks of 
the form A»B, with slightly adapted growth procedures. An advantage of this method would be 
a speedup by roughly a factor d in the large sparse eigensolver; for example, the application of h 
to a state in Eq. (fT2l) would then lead to 0{D^d) operations. In the infinite-system algorithm an 
obvious disadvantage would be that superblock lengths oscillate between odd and even; in the 
finite-system algorithm the question of the relative merits is much more interesting and will be 
discussed at length in Section l674l 

Obviously, for D — > oo, no truncations occur and DMRG becomes exact; increasing D re- 
duces truncation and therefore monotonically improves observables, which one extrapolates in 
D — > oo (even better, in the truncation error e — > 0, for which local observables often show 
effectively linear error dependence on e) for optimal results. 

For more details on DMRG and its applications, 1 refer to 10] • 

3. DMRG and entanglement: why DMRG works and why it fails 

The DMRG algorithm quite naturally leads to the consideration of bipartite quantum systems, 
where the parts are A» and "B. For an arbitrary bipartition, = 2^ 4'ij\i)A\i)B^ where the states 
|/)^ and form orthonormal bases of dimensions Na and A^^ respectively. Thinking of the i/r/^ 
as entries of a rectangular matrix *F (dimension A^^ x Nb), the reduced density matrices pa and 
Pb take the form 

PA = ^P*!'' PB = *1''*1'. (16) 

If we assume that we know exactly, but can approximate it in DMRG only with at most D 
states per block, the optimal DMRG approximation is provided by retaining as block states the 
eigenstates belonging to the D largest eigenvalues. If we happen to know the eigenspectra of 
reduced density operators of we can easily assess the quality a DMRG approximation can 
have; it simply depends on how quickly the eigenvalues Wa decrease. In fact, such analyses have 
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been carried out for some exactly solved systems in one and two dimensions II85L 1861 1871 18811891] . 
They reveal that in one dimension for gapped systems eigenvalues Wg generically decay expo- 
nentially fast (roughly as e which explains the success of DMRG, but in two-dimensional 
stripe geometries of size LxW, L'» W, c cc l/W, such that with increasing width W (increasing 
two-dimensionality) the eigenspectrum decay is so slow as to make DMRG inefficient. 

Usually, we have no clear idea about the eigenvalue spectrum; but it turns out that in such 
cases entanglement entropies can serve as "proxy" quantities, namely the von Neumann entan- 
glement or entanglement entropy. It is given by the non-vanishing part of the eigenvalue spectrum 
of Pa (identical to that of ps, as we will discuss below) as 



-TrpA log2 Pa = - ^ Wfl log2 w^. (17) 



It would seem as if we have gained nothing, as we don't know the Wg, but general laws about 
entanglement scaling are available. If we consider a bipartitioning A|B where AB is in the 
thermodynamic limit and A of size L®, with D the spatial dimension, the so-called area laws 
1 31, 9^ 91, 92, 93] predict that for ground states of short-ranged Hamiltonians with a gap to 



excitations entanglement entropy is not extensive, but proportional to the surface, i.e. S(A\B) ~ 
as opposed to thermal entropy. This implies S ~ est. in one dimension and S ~ Lin two 
dimensions. At criticality, a much richer structure emerges: in one dimension, S - ^ log2 L + k, 
where c and c are the (an)holonomic central charges from conformal field theory [^29^ 30,1; in two 
dimensions, bosonic systems seem to be insensitive to criticality (i.e. S oc L) i9li 19411 . whereas 
fermionic systems get a logarithmic correction S cc L log2 L for a one-dimensional Fermi surface 
(with a prefactor proportional to its size), but seem to grow only sublogarithmically if the Fermi 
surface consists of points ll94ll95ll . It should be emphasized that these properties of ground states 
are highly unusual: in the thermodynamic limit, a random state out of Hilbert space will indeed 
show extensive entanglement entropy with probability 1 . 

In a mathematically non-rigorous way one can now make contact between DMRG and the 
area laws of quantum entanglement: between two D-dimensional state spaces for A and B, the 
maximal entanglement is log2 D in the case where all eigenvalues of pA are identical and D^^ 
(such that Pa is maximally mixed); meaning that one needs a state of dimension 2^ and more 
to encode entanglement S properly. This implies that for gapped systems in one dimension an 
increase in system size will not lead to a strong increase in D; in two dimensions, D ~ 2^, 
such that DMRG will fail even for relatively small system sizes, as resources have to grow 
exponentially (this however does not exclude very precise results for small two-dimensional 
clusters or quite large stripes). Critical systems in one dimension are borderline cases: D ~ ; 
this means that the thermodynamic Umit is not reachable, but the growth of D is sufficiently 
slow (usually the power is weak, say 1/3 or 1/6, due to typical values for central charges) such 
that large system sizes (L ~ 0(10"')) can be reached; this allows for very precise finite-size 
extrapolations. 

Obviously, this argument implicitly makes the cavalier assumption that the eigenvalue spec- 
trum is close to flat, which leads to maximal entanglement, such that an approximate estimation 
of D can be made. In practice, the spectrum is dictated by the problem and indeed far from flat: 
as we have seen, it is in fact usually exponentially decaying. But numerically, it turns out that for 
standard problems the scaling of the resource D is predicted correctly on the qualitative level. 

It even turns out that in a mathematically strict analysis, von Neumann entanglement does not 
allow a general prediction of resource usage: this is because one can construct artificial eigen- 
value spectra that allow or forbid efficient simulation, while their von Neumann entanglement 
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would suggest the opposite, following the above argument ll96ll . Typical many-body states of in- 
terest, however, do not have such "pathological" spectra. In fact, Renyi entanglement entropies, a 
generalization of von Neumann entanglement entropies, do allow mathematically rigourous con- 
nections 1 96], but usually are hard to calculate, with criticality in one dimension as an exception 
due to conformal field theory. 

4. Matrix product states (MPS) 

If we consider our paradigmatic problem, the one-dimensional Heisenberg antiferromagnet, 
the key problem is that Hilbert space seems to be exponentially big (d^ = 2^). Looking for the 
ground state may therefore seem like looking for a needle in the haystack. The claim is that at 
least for local Hamiltonians with a gap between ground state and first excited state, the haystack 
is not very big, effectively infinitesimally small compared to the size of the full Hilbert space, as 
we have already seen from the very peculiar entanglement scaling properties. What is even more 
important, this relevant corner of Hilbert space can be parametrized efficiently, i.e. with modest 
numerical resources, operated upon efficiently, and efficient algorithms of the DMRG type to 
solve important questions of quantum physics do exist. This parametrization is provided by the 
matrix product states (MPS). 

Maybe the two DMRG algorithms explained above seem to be very cumbersome to imple- 
ment. But it turns out that if we do quantum mechanics in the restricted state class provided by 
matrix product states, DMRG and other methods almost force themselves on us. The manipula- 
tion of matrix product states seems to be very complicated at first, but in fact can be formalized 
beautifully, together with a graphical notation that allows to generate permitted operations almost 
automatically; as any good formalism (such as bra and ket), it essentially enforces correctness. 

4.1. Introduction of matrix product states 

4.1.1. Singular value decomposition (SVD) and Schmidt decomposition 

Throughout the rest of this paper, we will make extensive use of one of the most versatile 
tools of linear algebra, the so-called singular value decomposition (SVD), which is at the basis 
of a very compact representation of quantum states living in a bipartite universe AB, the Schmidt 
decomposition. Let us briefly recall what they are about. 

SVD guarantees for an arbitrary (rectangular) matrix M of dimensions (A^^ xNe) the existence 
of a decomposition 

M^USV\ (18) 

where 

• f/ is of dimension (A^^ x vnm{NA,NB)) and has orthonormal columns (the left singular 
vectors), i.e. U^U - /; if Na < Nb this implies that it is unitary, and also C/C/' = /. 

• 5 is of dimension {min(NA,NB) x mm{NA,NB)), diagonal with non-negative entries Saa = 
Sg. These are the so-called singular values. The number r of non-zero singular values is 
the (Schmidt) rank of M. In the following, we assume descending order: si > S2 > ■ . ■ > 

Sr > 0. 

• y ' is of dimension {mm{NA,NB) x A^^) and has orthonormal rows (the right singular vec- 
tors), i.e. y ' y - I. If Na > Nb this implies that it is unitary, and also yy^ = /. 
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Figure 4: Resulting matrix shapes from a singular value decomposition (SVD), coiresponding to the two rectangular 
shapes that can occur. The singular value diagonal serves as a reminder that in M = USV ^ S is purely non-negative 
diagonal. 



This is schematically shown in Fig.|4] Singular values and vectors have many highly interesting 
properties. One which is of practical importance in the following is the optimal approximation 
of M (rank r) by a matrix M' (with rank r' < r) in the Frobenius norm ||M||^ - ^,7 lA^/jP induced 
by the inner product {M\N) - Tr M'^'N. It is given by 



M'^US'V'^' S' ^diag(si,S2,. 



sO,...), 



(19) 



i.e. one sets all but the first r' singular values to be zero (and in numerical practice, will shrink 
the column dimension of U and the row dimension of V ' accordingly to r')- 

As a first application of the SVD, we use it to derive the Schmidt decomposition of a general 
quantum state. Any pure state \t//} on AB can be written as 



(20) 



where and {\j}g} are orthonormal bases of A and B with dimension A^^ and Nb respectively; 
we read the coefficients as entries of a matrix ^F. From this representation we can derive the 
reduced density operators pA - TrB|i/r)(i/'| and ps - TrA\>ff){^\, which expressed with respect to 
the block bases take the matrix form 



PA = pB = ^P'^F. 

If we carry out an SVD of matrix ^' in Eq. ( l20b . we obtain 

min(A',i,A'B) 

'■^^ = Z Tj UiaSaaV*jJi}A\j}B 
ij a=l 

fl=l V I V 
min(A',i,A'B) 



(21) 



^ Sa\a}A\a}B. 



(22) 



a=l 



Due to the orthonormality properties of U and the sets {la)^} and are orthonormal and 
can be extended to be orthonormal bases of A and B. If we restrict the sum to run only over the 
r < mm{NA, Nb) positive nonzero singular values, we obtain the Schmidt decomposition 



|iA) = ^ Sa\a)A\a)B. 



(23) 



£1=1 
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It is obvious that r = I corresponds to (classical) product states and r > 1 to entangled (quantum) 
states. 

The Schmidt decomposition allows to read off the reduced density operators for A and B 
introduced above very conveniently: carrying out the partial traces, one finds 

r r 

showing that they share the non-vanishing part of the spectrum, but not the eigenstates. The 
eigenvalues are the squares of the singular values, Wa - s^, the respective eigenvectors are the 
left and right singular vectors. The von Neumann entropy of entanglement can therefore be read 
off directly from the SVD, 

r 

Sa\b{W)) = -Tr Pa logzPA = - ^ loga s^. (25) 

a=l 

In view of the large size of Hilbert spaces, it is also natural to approximate li/r) by some 
spanned over state spaces of A and B that have dimension r' only. This problem can be related 
to SVD, because the 2-norm of is identical to the Frobenius norm of the matrix ^P, 

ij 

if and only if the sets {|/)) and are orthonormal (which is the case here). The optimal ap- 
proximation is therefore given in the 2-norm by the optimal approximation of by V in the 
Frobenius norm, where 'i' is a matrix of rank r'. As discussed above, V = US'V'' , where 
S' - diag(ii, . . . , irs 0, . . .), constructed from the largest singular values of *P. Therefore, the 
Schmidt decomposition of the approximate state reads 

where the Sa must be rescaled if normaUzation is desired. 

4.1.2. QR decomposition 

While SVD will be seen to cover all our needs, sometimes it is an overkill: in many cases of 
the expression M - USV\ we are only interested in the property U^U = I and the product S V^, 
for example whenever the actual value of the singular values will not be used explicitly. Then 
there is a numerically cheaper technique, QR decomposition, which for an arbitrary matrix M of 
dimension {Na x Nb) gives a decomposition 

M = QR, (28) 

hence the name, where Q is of dimension {Na X Na) and unitary, Q^Q = I = QQ\ and R is of 
dimension (Na X Nb) and upper triangular, i.e. Rij = if j > j. This full QR decomposition can 
be reduced to a thin QR decomposition: assume Na > Nb'. then the bottom Na - Nb rows of R 
are zero, and we can write 



M=Q 



Ri 
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Ri 




QiRu (29) 



where Qi is now of dimension (A^^ x Nb), Ri of dimension (Nb X Nb), and while Q\Qi - I, in 
general QiQ[ I. Whenever I will refer to a QR decomposition in the following, I will imply 
the thin one. It should also be clear that the matrices Q (or Qi) share properties with U from 
SVD, but are not the same in general; but as we will see that the MPS representation of states is 
not unique anyways, this does not matter. 

4.1.3. Decomposition of arbitrary quantum states into MPS 

Consider a lattice of L sites with i/-dimensional local state spaces {cr;) on sites i = 1, . . . ,L. 
In fact, while we will be naturally thinking of a one-dimensional lattice, the following also holds 
for a lattice of arbitrary dimension on which sites have been numbered; however, MPS generated 
from states on higher-dimensional lattices will not be manageable in numerical practice. The 
most general pure quantum state on the lattice reads 

ki,...,cri), (30) 

o-i,...,o-i 

where we have exponentially many coefficients Co-, with quite oblique content in typical quan- 
tum many-body problems. Let us assume that it is normalized. Can we find a notation that gives 
a more local notion of the state (while preserving the quantum non-locality of the state)? Indeed, 
SVD allows us to do just that. The result may look quite forbidding, but will be shown to relate 
profoundly to familiar concepts of quantum physics. There are three ways of doing this that are 
of relevance to us. 

( i) Left-canonical matrix product state. In a first step, we reshape the state vector with 
components into a matrix T of dimension (d X d^~^), where the coefficients are related as 

^£ri,(o"2---f^L) ~ ^0-i...(Ti- (31) 

An SVD of T gives 

ri ri 
Ctri...tri = *I'tri,(tr2...a-i) — Ucr^ai^ ai,aiiy^)aiXa-2-a-L) = f^cri,aiCaia-2...a-L> (32) 

where in the last equality S and V ' have been multiplied and the resulting matrix has been 
reshaped back into a vector. The rank is ri < d. We now decompose the matrix U into a 

collection of d row vectors A""' with entries A^' - U,r,,ai ■ At the same time, we reshape Cai(T2...iTL 
into a matrix ^(aio-2),(o-3- -tri) of dimension (rid x d^^^), to give 



c, 



(Ti...(Tl 



Z^ai'*^(ai(^2),(CT-3-o"i)- (33) 



^ is subjected to an SVD, and we have 

Ccri...crL — Z Z'^ai'^(aio-2),a2'^a2,a2(^^)a2,(a"3...cri) = ^al^al.a^ (aicri),(a-i-a-L)' (34) 

a\ ai fl] ^2 

where we have replaced t/ by a set of d matrices A""^ Qf dimension (ri x r2) with entries A'^l^^^ = 
Uiaia-2),a2 multiplied S and V^, to be reshaped into a matrix of dimension {rid x d^~^), 
where < r\d < Upon further SVDs, we obtain 

ai,...,aL-i 
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or more compactly 

= A'^'A-^^ . . . A'^^-'A'^S (36) 

where we have recognized the sums over oi, 02 and so forth as matrix multiplications. The last 
set of matrices A°'^ in fact consists of column vectors. If we wish, dummy indices 1 may be 
introduced in the first and last A to turn them into matrices, too. In any case, the (arbitrary) 
quantum state is now represented exactly in the form of a matrix product state: 

\ilf)= Yj A''^A''K..A''^'^A''^\ai,...,(TL). (37) 

cr\,...,(TL 

Let us study the properties of the A-matrices. The maximal dimensions of the matrices are 
reached when for each SVD done the number of non-zero singular values is equal to the upper 
bound (the lesser of the dimensions of the matrix to be decomposed). Counting reveals that the 
dimensions may maximally be (1 x J), {dX(f), (J^/^-i xd^'^), (d^'^xd^'^'^), (d^xd), {dx 
1), going from the first to the last site (1 have assumed L even for simplicity here). This shows 
that in practical calculations it will usually be impossible to carry out this exact decomposition 
expUcitly, as the matrix dimensions blow up exponentially. 

But there is more to it. Because at each SVD U^U = I holds, the replacement of f/ by a set 
of A"' entails the following relationship: 

at-i<Tt 

or, more succinctly, 

Y^A'^''^A"'' = I. (38) 

(Tf 

Matrices that obey this condition we will refer to as left-normalized, matrix product states that 
consist only of left-normalized matrices we will call left-canonical. In fact, a closer look reveals 
that on the last site the condition may be technically violated, but as we will see once we calculate 
norms of MPS this corresponds to the original state not being normaUzed to 1 . Let us ignore this 
subtlety for the moment. 

In view of the DMRG decomposition of the universe into blocks A and B it is instructive to 
spUt the lattice into parts A and B, where A comprises sites 1 through { and B sites {-t-l through 
L. We may then introduce states 

Ma = Yj ('^"''^"'•••'^"Oi.«,|cri,...,cr,) (39) 

o-i,...,a-, 

Mb = 2 (A'''*'A-'''---A'"-)a,.i\(rM,...,(rL) (40) 

a"f+i,-,o"L 

such that the MPS can be written as 

\0 = 2 l«f>Al«^>i*- (41) 

at 
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This pairing of states looks tantalizingly close to a Schmidt decomposition of \\p), but this is 
not the case. The reason for this is that while the {|a^)A) form an orthonormal set, the {Iflf)^) in 
general do not. This is an immediate consequence of the left-normality of the A-matrices. For 
part A we find 

A{a'i\ai)^ = ^ (A' 

o-i,...,o-, 

0-i,...,(T, 

o-i,...,o-f 

where we have iteratively carried out the sums over a-\ through crt and used left-normality. On 
the other hand, the same calculation for part B yields 

= 2 (A-^^-f . . . A'^'-'-'-)i,,-(A'^^^' . . . A-^O^,,! 

which cannot be simplified further because in general 2o- '4'^A'^' + I. 

The change of representation of the state coefficients can also be represented graphically 
(Fig- 121). Let us represent the coefficient Co-,...cri as a black box (with rounded edges), where the 
physical indices cri through cri stick out vertically. The result after the first decomposition we 
represent as in the second line, where we have on the left hand site an object representing A^,' , 
on the right Ca^a■2...a■L■ The auxiliary degrees of freedom (ai) are represented by horizontal lines, 
and the rule is that connected lines are summed over. The second step is then obvious, we have 
Afl|' , then A^^^, and on the right Ca2a-i...a-L^ with all connected lines summed over. In the end, we 
have arrived at L A-matrices multiplied together and labelled by physical indices (last line of the 
figure). 

The graphical rules for the A-matrices, that on the first and last site are row and column 
vectors respectively, are summarized in Fig. |6] a site t is represented by a solid circle, the 
physical index o-( by a vertical line and the two matrix indices by horizontal lines. 

Let me conclude this exposition by showing the generation of a left-canonical matrix product 
state by a sequence of QR decompositions. We start as 

C(Ti...o-i — *Po-i,(o"2...(Tl) — Qa-i,a\Rai,(a-2...crL) — ^l.m ^ (aiO'2),(o'y-crL)^ (42) 
fli fli 

where we reshape Q ^ A and 7? ^ *P in analogy to the SVD procedure. The next QR decompo- 
sition yields 

Ccri...tTL — ^ 6(flio-2),fl2''^fl2,(a"3--0"l) ~ ^\]a^^ai,a^(aiO'~MP'i--0'L) (43) 
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Figure 5: Graphical representation of an iterative construction of an exact MPS representation of an arbitrary quantum 
state by a sequence of singular value decompositions. 




Figure 6: Graphical representation of A-matrices at the ends and in the bulk of chains: the left diagram represents A 'j^J^^ , 
the row vector at the left end, the right diagram represents A'^j-^, the column vector at the right end. In the center there 
is 4 
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and so on (on the right half of the chain, thin QR is needed, as an analysis of the dimensions 
shows). Q' Q - I implies the desired left-normalization of the A-matrices. If numerically fea- 
sible, this is faster than SVD. What we lose is that we do not see the spectrum of the singular 
values; unless we use more advanced rank-revealing QR decompositions, we are also not able to 
determine the ranks r\,r2, ■ ■ ., unlike in SVD. This means that this decomposition fully exploits 
the maximal A-matrix dimensions. 

( ii) Right-canonical matrix product state. Obviously, there was nothing specific in the de- 
composition starting from the left, i.e. site 1. Similarly, we can start from the right in order to 
obtain 

Ca-i...a-L — *J'(o-i...o-L_i),o-i 

- ^ t^(a-i...o-i_i),ai_i'5fli_i,fli_i(^')at_i,tri 

~ ^(cr\...a-L-2),(o'L-iaL-\)^ai_i 

aL-\ 

~ ^(o-i...a-L-2),aL-i^ aL-2,aL-2(^'^)aL-2,{o-L-iaL-l)^ai-i 

aL-2,aL-\ 

~ ^ ^(o-\-a-L-Ma-L-2aL-2)^aiJ,ai_,^at-i ~ ■ ■ ■ 
aL-2,aL-\ 

/ 1 ai a[,a2 ai_2AL-l Ol-i 

ai,...,aL_i 

Here, we have reshaped (V^)ai_,.a-L into d column vectors -So/:,, (V^ )(fli_,o-i_i),fli_i into d ma- 
trices B^^'^o'.OL-i' ^^'^ ^° multiplied U and S before reshaping into *P at each step. 
The obvious graphical representation is given in Fig. [T] We do not distinguish in the graphical 
representation between the A- and Z?-matrices to keep notation simple. 

We obtain an MPS of the form 

^ B'^'B'^^..B'^^-'B'^^|c^l,...,c^i) (44) 

iTi,...,a-£, 

where the B-matrices can be shown to have the same matrix dimension bounds as the A matrices 
and also, from V V = /, to obey 

^B'^'B'^fl- = /, (45) 

such that we refer to them as right-normalized matrices. An MPS entirely built from such matri- 
ces we call right-canonical. 

Again, we can split the lattice into parts A and B, sites 1 through t and t + I through L, and 
introduce states 

\a(}A = Yj ('B"''B"'---'B"Oi,«,ki,...,crf) (46) 

o-i,...,crf 

\a()B = 2 {B^'-'B'''^^...B^^)a„i\crM,...,crL) (47) 

o-f+i,...,a-i 

such that the MPS can be written as 

m^Yj^ae)A\af)B. (48) 

Of 
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This pairing of states looks again tantalizingly close to a Schmidt decomposition of \tli), but this 
is again not the case. The reason for this is that while this time the {\ac)g] form an orthonormal 
set, the (Iflf)^) in general do not, as can be shown from the right-normality of the Z?-matrices. 

Again, the right-normalized form can be obtained by a sequence of QR decompositions. The 
difference to the left-normalized form is that we do not QR-decompose = QR, but T ' = QR, 
such that *P = R^Q^- This leads directly to the right-normalization properties of the Z?-matrices, 
if we form them from Q ' . Let me make the first two steps explicit; we start from 

Co-i.-a-L - *I'(o-i...o-i_,),o-i - ^(o-i ...o-,_| ),fl,_| QIl-,,o-l ~ ^ (o-i-o-L-2),(o-L-iaL-i)^ a^_,,\^ (49) 

reshaping ' into ^', 2 ' into B, and continue by a QR decomposition of *P' as 

Co-i...a-L - ^ ^lo-i...o-L-2),aL-2QaL-2,(o-L-iaL-l)^aL-i,l ~ ■•■0"i-3).(o-L-2flL-2)^fli J.Oi-l ^flt-l ,1 ' 

aL-\,aL-2 aL-\,aL-2 

(50) 

We have now obtained various different exact representations of in the MPS form, which 
already indicates that the MPS representation of a state is not unique, a fact that we are going to 
exploit later on. 

( Hi) Mixed-canonical matrix product state. We can also mix the decomposition of the state 
from the left and from the right. Let us assume we did a decomposition from the left up to site i, 
such that 

Co-i...o-f, — '...A ')a[S a(,a((y'^)a[,(a-c+\...crL)- (51) 

We reshape as *P(flfO-,+,...o-i_,),o-i and carry out successive SVDs from the right as in the original 
decomposition from the right, up to and including site crt+i, in the last SVD t/(afO-,+i),af+,'5af^,,flf^, 
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Figure 8: Graphical representation of an exact MPS obtained by a sequence of singular' value decompositions, start- 
ing from the left and right. The diamond represents the diagonal singular value matrix. Matrices to the left are left- 
normalized, to the right are right-normalized. 



remains, which we reshape to Bo/flf^i ■ Then we obtain 

flf+i,...,ai-i 

All Z?-matrices are right-normalized. This is simply due to the SVD for sites £ + 2 through L; on 
site { + 1, it follows from the property V ' V = /: 

^a[,a'^ — ^ (V )£if,(o-f+i...o-i)V(crf^,...o-i),£i^ 
o-f+i,... 

o-f+i,... 

where we use in the last line the right-normaUzation property of all the Z?-matrices on sites 
{ + 2, . . . , Lto obtain the desired result. 

We therefore end up with a decomposition 

= A'^' . . .A°''SB'^'^' . . . B'^S (53) 

which contains the singular values on the bond ^ -i- 1) and can be graphically represented as in 
Fig.i 

What is more, the Schmidt decomposition into A and B, where A runs from sites 1 to i and 
B from sites £ + 1 to L, can now be read off immediately. If we introduce vectors 

Ma = 2 (A- ...A'-Oi,«.ki,...,(rf) (54) 

o-i,...,o-f 

\ae}B = 2 ('B"'^'---'B"%,ik^+i'---'^L) (55) 

then the state takes the form (sa - Sa,a) 

\^) ^^Sa\ae)A\ae)B, (56) 

which is the Schmidt decomposition provided the states on A and are orthonormal respectively. 
But this is indeed the case by construction. 

(iv) Gauge degrees of freedom. By now, we have three different ways of writing an arbitrary 
quantum state as an MPS, all of which present advantages and disadvantages. While these three 
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are arguably the most important ways of writing an MPS, it is important to realise that the degree 
of non-uniqueness is much higher: MPS are not unique in the sense that a gauge degree of 
freedom exists. Consider two adjacent sets of matrices M'^< and M°"'+' of shared column/row 
dimension D. Then the MPS is invariant for any invertible matrix X of dimension (D x D) under 

M'^'^M"''X, M'^'*' M'^'*' . (57) 

This gauge degree of freedom can be used to simplify manipulations drastically, our three con- 
structions are just special cases of that. 

Several questions arise. Is there a connection between this notation and more familiar con- 
cepts from many-body physics? Indeed, there is a profound connection to iterative decimation 
procedures as they occur in renormalization group schemes, which we will discuss in Section 

Eiai 

The matrices can potentially be exponentially large and we will have to bound their size on 
a computer to some D. Is this possible without becoming too inaccurate in the description of the 
state? Indeed this is possible in one dimension: if we consider the mixed-canonical representa- 
tion, we see that for the exponentially decaying eigenvalue spectra of reduced density operators 
(hence exponentially decaying singular values Sg) it is possible to cut the spectrum following 
Eq. (IZTT i at the D largest singular values (in the sense of an optimal approximation in the 2-norm) 
without appreciable loss of precision. This argument can be generalized from the approximation 
incurred by a single truncation to that incurred by L - 1 truncations, one at each bond, to reveal 
that the error is at worst Wh 

L 

III<A)-I<^trunc)li <22e,(£'), (58) 

where ei{D) is the truncation error (sum of discarded squared singular values) at bond / incurred 
by truncating down to the leading D singular values. So the problem of approximability is as in 
DMRG related to the eigenvalue spectra of reduced density operators, indicating failure in two 
dimensions, and (a bit more tenuously) to the existence of area laws. 

4.1.4. MPS and single-site decimation in one dimension 

In order to connect MPS to more conventional concepts, let us imagine that we set up an 
iterative growth procedure for our spin chain, ^ — > ^ -H 1, as illustrated in Fig. |9] such that 
associated state spaces grow by a factor of d at each step. In order to avoid exponential growth, 
we now demand that state space dimensions have a ceiling of D. Once the state space dimension 
grows above D, the state space has to be truncated down by some as of now undefined ^rocsdms. 

Assume that somehow we have arrived at such a D-dimensional effective basis for our system 
(or left block A, in DMRG language) of length {-I, {|flf-i)A)- If the D basis states of the (left) 
block A of length { after truncation are {laf)^) and the local states of the added site {|crf)), we 
must have 

\at)A = ^ A{at-icre\at)A \ac-\)A\'^t) (59) 

flf-io-f 

for these states, with A{ai-icri\ac)A as of now unspecified. We now introduce at site { d matrices 
^[f]o-f of dimension {D x D) each, one for each possible local state |crf). We can then rewrite 
Eq. dig as 

\at)A = Yu ^^Z.\''t-^)A\(rt) (60) 

flf_l(Tf 
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\af)A 



1 £-11 1 i 

Figure 9: A block of length f - 1 is grown towards the right to a block of length ( by adding a site ( 





Figure 10: Graphical representation of A-matrices: the left diagram represents A a[_i,a[' 'he right diagram the conjugate 
A IJj'* . The solid circle represents the lattice sites, the vertical line the physical index, the horizontal lines the matrix 
indices. 



where the elements of the matrices A'^^'^' are given by (see Fig.fTOll 

= A{a(-icri\ai}^. (61) 

Let us make a short remark on notations right here: in A^^^°'' , [{] indicates which set of A-matrices 
is considered, and cr( which A-matrix in particular In the present case, this is a notational 
overkill, because the local state \o-(} is taken from the site where the matrices were introduced. 
In such cases, we drop one of the two {, usually [£]: 

We will, however, encounter situations where matrices A are selected by local states not on the 
site where they were introduced. In such cases, the full notation obviously has to be restored! 

Similarly, we will shorten |flf)^ — > \a(), when the fact that the state lives on block A is 
irrelevant or totally obvious. 

The advantage of the matrix notation, which contains the decimation procedure yet unspeci- 
fied, is that it allows for a simple recursion from a block of length f to the smallest, i.e. vanishing 
block. Quantum states obtained in this way take a very special form: 

- Z Z A- A-,...A- ,„,k,)k2)...|cr,) 

= J]iA-^'A-'\..A''')i,a,\cri}\<T2}...\ae}, (63) 

where i runs through all sites of block A. The index 'A' indicates that we are considering states on 
the left side of the chain we are building. On the first site, we have, in order to keep in line with 
the matrix notation, introduced a dummy row index 1 : the states of block length 1 are built from 
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Figure 11: Graphical representation of the recursive construction of a state \a() by contraction (multiplication) of A- 
matrices. Contractions run over all connected legs. 



I«£>i 



L 



i+l i+2 



Figure 12: A block B of length L - f - 1 is grown towards the left to a block B of length L — Chy adding site C +\. 

the local states on site 1 and the block states of the "block" of length 0, for which we introduce 
a dummy state and index 1. This means that A""' is in fact a (row) vector (cf. Fig.|6l). We also 
see that the left or row index of A correspond to states "left" of those labelled by the right or 
column index. Quite generally, we can show this construction as in Fig.[TT] if - as before - we 
introduce the rule that all connected legs are summed over (contracted). The advantage of the 
matrix notation is that we can hide summations in matrix multiplications. 

Similarly, we can build blocks to grow towards the left instead of to the right (Fig.fTSI): we 
have 

= ^ B<flf+lO-f+i|flf)B |flf+i)B|cr;+l) (64) 



or 



\ai)B- Yj <:«'iri«f+i)Bi^f+i) 



(65) 



with 



(66) 



We call matrices B to indicate that they emerge from a growth process towards the left, in DMRG 
language this would mean block B. Recursion gives 



\ae)B = ^(B'^^'B"'^^ . . . B"0o,,„ikf+i)|crf+2) ■ • ■ \<Tl), 



(67) 



where / runs from { + \ io L, the sites of block B. A similar dummy index as for position 1 is 
introduced for position L, where the B-matrix is a (column) vector. 

Note a slight asymmetry in the notation compared to the A-matrices: in order to be able to 
match blocks A and B later, we label block states according to the bond at which they terminate: 
bond £ connects sites £ and i+l, hence a labeling as in Fig.[T3] 

If we introduce A-matrices and B-matrices in this way, they can be seen to have very special 
properties. If we consider the growth from the left, i.e. A-matrices, and demand reasonably 
that the chosen states should for each block length be orthonormal to each other, we have using 
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Figure 13: Blocks A (sites 1 through t) and B (sites ( + I through L) are joined at bond f. States are labelled \a f)^ and 




results. 



Eq. ^ 

cr'i,cri a[ ^,ar-i 

= I^I^K'XAZ.. - Y^A-''A'^%,.„,. (69) 

Summarizing we find that the A-matrices are left-normalized: 

^A'^'A'^ = /. (70) 

A graphical representation is provided in Fig. [141 The multiplication can also be interpreted as 
the contraction of A and A* over both cr and their left index. 

Similarly, we can derive for B-matrices of blocks B built from the right that the right- 
normalization identity 

(71) 

holds (usually, A and B will be used to distinguish the two cases). See Fig. [15] This means that 
orthonormal states can always be decomposed into left- or right-normalized matrices in the MPS 
sense and that all states constructed from left- or right-normalized matrices form orthonormal 
sets, provided the type of normalization and the direction of the growth match. 

Let us take a closer look at the matrix dimensions. Growing from the left, matrix dimensions 
go as (1 X d), id X d^), (d^ x d^), (d^ x D), where I have assumed that d'^ > D. Then they continue 
at dimensions (D x D). At the right end, they will have dimensions (D x D), (D x d^), (d^ x d^), 
(d^ X d) and (d x 1). 

We can now again write down a matrix product state. Putting together a chain of length L 
from a (left) block A of length i (sites 1 to {) and a (right) block B of length L - ( (sites ^ -H 1 to 
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Figure 15: If two rigfti-normalized B-matrices are contracted over their right index and the physical indices, a 6 a'^ aj line 
results. 




Figure 16: Representation of an open boundary condition MPS. 



L), we can form a general superposition 

Inserting the states explicitly, we find 

1^^) = }](A'^' . . .A'^')i,«,^«„a;(B""' . . (73) 
tr 

The bold-faced cr stands for all local state indices, \cr) - \o-i,cr2, . . . , cr^). The notation suggests 
to interpret as a matrix; then the notation simplifies to 

= ^ A"^' . . .A'''^'B'^'*' . . .B'^^lcr). (74) 
cr 

If we allow general matrices and don't worry about left, right or no normalization, we can 
simply multiply the ^P-matrix into one of the adjacent A or B matrices, such that the general MPS 
for open boundary conditions appears (see Fig.fTSli: 

I,/,) = ^ M"^' . . . M'^' |cr) (MPS for OBC), (75) 
cr 

where no assumption about the normalization is implied (which is why I call matrices M). Due 
to the vectorial nature of the first and last matrices the product results in a scalar. This is exactly 
the form of an MPS already discussed in the last section. 

At this point it is easy to see how a matrix product state can exploit good quantum numbers. 
Let us focus on magnetization and assume that the global state has magnetization M. This 
Abelian quantum number is additive, M - XiM- We choose local bases {o",) whose states are 
eigenstates of local magnetization. Consider now the growth process from the left. If we choose 
the states to be eigenstates of local magnetization (e.g. by taking just the |cri)), then Eq. (|59] l 

29 




Figure 17: Representation of an open boundary condition MPS with good (additive) quantum numbers. Physical states 
and bonds become directed, such that the quantum numbers on the ingoing lines equal those on the outgoing lines. For 
the dummy bonds before the first and after the last site we set suitable values to fix the global good quantum number. 



allows us to construct by induction states \ac) that are eigenstates of magnetization, provided the 
matrices Aa'^_^ af obtain a block structure such that for each non-zero matrix element 

M(\ae^i)) + MQcre)) ^ M(\ae)) (76) 

holds. This can be represented graphically easily by giving directions to the lines of the graphical 
representation (Fig.fTTli. with ingoing and outgoing arrows. The rule is then simply that the sum 
of the magnetizations on the ingoing lines equals that on the outgoing lines. In order to enforce 
some global magnetization M, we may simply give magnetization values and M to the ingoing 
and outgoing dummy bonds before the first and after the last site. We may envisage that the 
indices of the MPS matrices are multiindices for a given magnetization allowing degeneracy, 
leading to elegant coding representation. An inversion of the bond arrows would directly tie in 
with the structure of B-matrices from the growth from the right, but proper book-keeping gives 
us lots of freedom for the arrows: an inversion means that the sign has to be reversed. 

In order to use good quantum numbers in practice, they have to survive under the typical 
operations we carry out on matrix product states. It turns out that all operations that are not 
obviously unproblematic and maintain good quantum numbers can be expressed by SVDs. An 
SVD will be applied to matrices like A(ai_,a-,),ar If we group states |fl,_icr,) and |fl,) according to 
their good quantum number, A will consist of blocks; if we rearrange labels appropriately, we 
can write A = Ai®A2®. . . - UiSiVl@U2S2Vl®. . . orA = USV^' where U - t/i et/2©. - . and 
so forth. But this means that the new states generated from |fl,_icr,) via U will also have good 
quantum numbers. When the need for truncation arises, this property of course still holds for 
the retained states. If we replace SVD by QR where possible and carry it out on the individual 
blocks. A, = QiRi, the unitary matrices Qi transform within sets of states of the same quantum 
numbers, hence they remain good quantum numbers. 

Let us now assume that our lattice obeys periodic boundary conditions. At the level of the 
state coefficients Co-, o-^ there is no notion of the boundary conditions, hence our standard form 
of an MPS is capable to describe a state that reflects periodic boundary conditions. In that 
sense it is in fact wrong to say that Eq. (|75] | holds only for open boundary conditions. It is 
true in the sense that the anomalous structure of the matrices on the first and last sites is not 
convenient for periodic boundary conditions; indeed, the entanglement across the bond (L, 1) 
must be encoded as stretching through the entire chain. This leads to numerically very inefficient 
MPS representations. 

For periodic boundary conditions the natural generalization of the MPS form is to make aU 
matrices of equal dimensions (D x D); as site L connects back to site 1, we make the MPS 
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Figure 18: Representation of a periodic boundary condition MPS; the long line at the bottom coiresponds to the trace 
operation. 

consistent with matrix multiplications on all bonds by taking the trace (see Fig.fTSli: 

= Tr(M'^' ■ ■ ■ M'^^)\(r) (MPS for PBC). (77) 

While a priori not more accurate than the other, it is much better suited and computationally far 
more efficient. 

In this section, our emphasis has been on approximate representations of quantum states 
rather than on usually unachievable exact representations. While we have no prescription yet 
how to construct these approximate representations, some remarks are in order. 

Even an approximate MPS is still a linear combination of all states of the Hilbert space, 
no product basis state has been discarded. The limiting constraint is rather on the form of the 
linear combinations: instead of coefficients, dL matrices of dimension (D x D) with a matrix- 
valued normalization constraint that gives LD^ scalar constraints have (d - l)LD^ independent 
parameters only, generating interdependencies of the coefficients of the state. 

The quality of the optimal approximation of any quantum state for given matrix dimensions 
will improve monotonically with D: take D < D', then the best approximation possible for D 
can be written as an MPS with D' with {D x D) submatrices in the (D' x D') matrices and all 
additional rows and columns zero. They give further parameters for improvement of the state 
approximation. 

Product states (with Schmidt rank 1 for any Schmidt decomposition) can be written exactly 
using D - I MPS. Real quantum physics with entangled states starts at D - 2. Given the 
exponential number of coefficients in a quantum state, it may be a surprise that even in this sim- 
plest non-trivial case interesting quantum physics can be done exactly! But there are important 
quantum states that find a compact exact expression in this new format. 

4.1.5. The AKLT state as a matrix product state 

In order to make the MPS framework less abstract, let us construct the MPS representation of 
a non-trivial quantum state. One of the most interesting quantum states in correlation physics is 
the Affleck-Kennedy-Lieb-Tasaki state introduced in 1987, which is the ground state of the AKLT 
HamiltonianiUlMl 

^ = J]S,-S,+, + i(S,-S,-+i)2, (78) 

where we deal, exceptionally in this paper, with S - I spins. It can be shown that the ground state 
of this Hamiltonian can be constructed as shown in Fig. [19] Each individual spin-1 is replaced by 
a pair of spin-j which are completely symmetrized, i.e. of the four possible states we consider 
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spin-1 singlet spin- 1/2 

Figure 19: The Affleck-Kennedy-Lieb-Tasaki (AKLT) state is built from expressing the local spin-1 as two totally sym- 
metrized spin- J particles which are linked across sites by singlet states. The picture shows the case of PBC, for OBC the 
"long" bond is cut, and two single spins- ^ appear at the ends. 



only the three triplet states naturally identified as 5 = 1 states: 

1+) = ITT) 

I Ti> + I iT> 



10) = 



V2 



(79) 



I-) = lU) 

On neighbouring sites, adjacent pairs of spin-j are linked in a singlet state 

I Ti) - I iT) 
V2 ■ 



(80) 



As it turns out, this state can be encoded by a matrix product state of the lowest non-trivial 
dimension D = 2 and contains already lots of exciting physics ll33ll34i 13511 . In the language of 
the auxiliary 2L spin-| states on a chain of length L any state is given as 



X '■"blab) 

a b 



(81) 



with |a) = |fli , . . . , ai) and |b) = \bi, . . ., bi) representing the first and second spin-| on each 
site. We now encode the singlet bond S, on bond / connecting sites / and / H- 1 as 

|2W) = 2 2/,J^-)|fl,+i) (82) 



introducing a 2 x 2 matrix E 



in 



Then the state with singlets on all bonds reads 



IfAi;) = ^b^a.'^bia, ■ ■ ■ l^ht-xaj^btai |ab) 

a b 



(83) 



(84) 



for periodic boundary conditions. If we consider open boundary conditions, S/j^oj is omitted and 
the first and last spin-:^ remain single. 

Note that this state is a product state factorizing upon splitting any site / into its two con- 
stituents. We now encode the identification of the symmetrized states of the auxiliary spins 
with the physical spin by introducing a mapping from the states of the two auxiliary spins- i, 
\ai)\bi) e {| T)J i))®^ to the states of the physical spin-1, |cr,) € |0), |-)). To represent 
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Eq. ( |79] |. we introduce M'^^i^\a-){ab\, with \ab) and \cr) representing the auxihary spins and the 
physical spin on site /. Writing M^"^^ as three 2x2 matrices, one for each value of |cr), with rows 
and column indices standing for the values of \a) and \b), we find 



10 
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V2 
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V V2 





1 



The mapping on the spin-1 chain Hilbert space {\cr)) then reads 



y y M'^' M'\ . ..M'^', |o->(ab|. 

/ I / I IIlOi 0202 flLOL 



(T ab 



li^s) therefore is mapped to 



(T ab 



or 



|(A> = Tr(M'^'SM°'=S . . . M°'^S])|o-), 



using the matrix notation. To simplify further, we introduce A°" = M°"S, such that 



i 







2 

+^ 





The AKLT state now takes the form 

|^) = ^Tr(A'^■A"^..A'^^)|o■). 







(85) 

(86) 

(87) 
(88) 

(89) 
(90) 



cr 



Let us left-normalize the A"". 2o-A°"'^A°" - |/, which implies that the matrices A should be 
rescaled by such that we obtain normalized matrices A, 



A+ = 


Vi 


A» = 


I 




1 


A- = 




l-Vf 















VI J 










(91) 



This normalizes the state in the thermodynamic limit: we have 

(^1^) = ^ Tr(A'^' . . . A'^0*Tr(A'^' . . . A'^O 



cr 

= Tr 



^A°"'*®A'"' ... ^A'^'*®A'^^ 



In this expression, the /I, ai^e the 4 eigenvalues of 



= 2- 



i 



-J 



1 
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4 



(92) 




111") 



Figure 20: Overlap between two states \(fi) and All contractions (sums) over same indices are indicated by arrows. 



namely But then <(A|(A> = 1 + 3(-l/3)^ ^ 1 fori ^ oo. 

The methods of the next section can now be used to work out analytically the correlators of 
the AKLT state: antiferromagnetic correlators are decaying exponentially, {S^-S'j} oc (-1/3)'"-', 

whereas the string correlator {S^e"'^''''''^'i:S^j} = -4/9 for j - i > 2, indicating hidden order 

To summarize, it has been possible to express the AKLT state as a D = 2 matrix product 
state, the simplest non-trivial MPS! In fact, the projection from the larger state space of the 
auxiliary spins which are linked by maximally entangled states (here: singlets) onto the smaller 
physical state space can also be made the starting point for the introduction of MPS and their 
higher-dimensional generalizations ifsH 671. 



4.2. Overlaps, expectation values and matrix elements 

Let us now turn to operations with MPS, beginning with the calculation of expectation values. 
Expectation values are obviously special cases of general matrix elements, where states and 
10) are identical. Staying in the general case, let us consider an overlap between states li/r) and 
10), described by matrices M and M, and focus on open boundary conditions. 

Taking the adjoint of |0), and considering that the wave function coefficients are scalars, the 
overlap reads 

(0|i/,) = ^ M'"' * . . . M'^^*M'^' ...M'^'-. (93) 
cr 

Transposing the scalar formed from the M . . .M (which is the identity operation), we arrive at 
adjoints with reversed ordering: 

(0|(/,) = ^ M'^it . . . M'^' "'m'"' ...M'^^. (94) 
cr 

In a pictorial representation (Fig. l20b . this calculation becomes much simpler, if we follow the 
rule that all bond indices are summed over. 



4.2.1. Efficient evaluation of contractions 

Evaluating expression ( l94b in detail shows the importance of finding the right (optimal) order 
of contractions in matrix or more generally tensor networks. We have contractions over the 
matrix indices implicit in the matrix multiplications, and over the physical indices. If we decided 
to contract first the matrix indices and then the physical indices, we would have to sum over 
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Figure 21: Overlap between two states \(fi) and with indication of the optimal sequence of contractions, running like 
a zipper through the chain. 



Strings of matrix multiplications, which is exponentially expensive. But we may regroup the 
sums as follows: 



(01,/,) = ^ M'^it ... Yj ^"'^ 



(95) 



This means, in the (special) first step we multiply the column and row vectors M°"' ' and M°"' to 
form a matrix and sum over the (first) physical index. In the next step, we contract a three-matrix 
multiplication over the second physical index, and so forth (Fig.[2TJ. The important observation 
is that from the second step onwards the complexity does not grow anymore. Also, it is of 
course most efficient to decompose matrix multiplications ABC as {AB)C or A(BC). Then we 
are carrying out (2L - \)d multiplications, each of which is of complexity O(D^), ignoring for 
simplicity that matrices are non-square in general at the moment. The decisive point is that we 
go from exponential to weak polynomial complexity, with total operation count O(LD^d). 

What is also immediately obvious, is that for a norm calculation and OBC having a state 
in left- or right-normalized form immediately implies that it has norm 1 . In the calculation above 
it can be seen that for left-normalized matrices A, the innermost sum is just the left-normahzation 
condition, yielding /, so it drops out, and the next left-normalization condition shows up, until 
we are through the chain (Fig. 



\ a-i J 
\ \ 



^A'"^"'" ... ^A'"^"'" 2 A'"' "'"A 
^A'"^"'" ... ^A'^^'A' 



A"^ 



To calculate general matrix elements, we consider {(p\(P^G^'^ . . . tensored operators act- 
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Figure 22: Steps of a norm calculation for a left normalized state \if/) by subsequent applications of the contraction rule 
for left-normalized A-matrices. 




Figure 23: Matrix elements between two states \(p) and \tfr) are calculated like the overlap, with the operators inserted at 
the right places, generating a double sum of physical indices there, as indicated by the arrows. 



ing on sites / and j. The matrix elements of such operators are taken from 



(96) 



Let us extend this to an operator on every site, which in practice will be the identity on almost 
all sites, e.g. for local expectation values or two-site correlators. We are therefore considering 
operator matrix elements O'^' 0°'-'°^^- ■ ■ ■ 0"''-'°''l . In the analytical expression, we again transpose 
and distribute the (now double) sum over local states (matrix multiplications for the M-matrices 
are as before): 

^ M""' * . . . M°"'^*0°'' 0°'--'^^- ■ ■ ■ 0'^''-°'lM°''' . . . M'^L 

(T,(T' 



Vo-j ,0-, 















j 







This amounts to the same calculation as for the overlap, with the exception that formally the 
single sum over the physical index turns into a double sum (Fig.l23ll. For typical correlators the 
double sum will trivially reduce to a single sum on most sites, as for most sites only the identity 
acts, O'^'^ = /; on the few non-trivial sites, of the up to S matrix elements, most will be zero for 
conventional operators, strongly restricting the number of terms, so essentially the operational 
count is 0(LD^d) again. 

Important simplifications for expectation values {t//\0^^^\ifr} are feasible and should be ex- 
ploited whenever possible: Assume that we look at a local operator O^'^^ and that normalizations 
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Figure 24: for a state with left- and riglit-normalized matrices to tiie left and riglit of site £. 



are such that to the left of site £ all matrices are left-normalized and to the right of site £ all matri- 
ces are right-normalized; the status of site £ itself is arbitrary. Then left- and right-normalization 
can be used to contract the network as in Fig. l22l without explicit calculation, such that just two 
matrices remain (Fig.l24li. The remaining calculation is just 

((/r|Ot^V> = X C'^"^'Tr(M'^'"'"M'"0, (97) 

an operation of order O(D^d^), saving one order of D in calculation time. As we will encounter 
algorithms where the state is in such mixed-canonical representation, it makes sense to calculate 
observables "on the sweep". This is just identical to expectation values on the explicit sites of 
DMRG. 

A frequently used notation for the calculation of overlaps, expectation values and matrix 
elements is provided by reading the hierarchy of brackets as an iterative update of a matrix, 
which eventually gives the scalar result. We now introduce matrices C'^', where C'"' is a dummy 
matrix, being the scalar 1. Then an overlap {(pli//) can be carried out iteratively by running { from 
1 through L: 

Cl^l = ^M'^f"'"C^^"'lM'^^ (98) 

where C^^^ will be a scalar again, containing the result. For operators taken between the two 
states, the natural extension of this approach is 

C'^J = ^ 0'"''^''M'^'^C^^-^^M<. (99) 

0"f,0"'f 

Again, the right order of evaluating the matrix products makes a huge difference in efficiency: 



(100) 



reduces an operation 0{D'^(f ) to OiD^d) + 0{D^cf ) + O(D^d). 

Of course, we can also proceed from the right end, introducing matrices D^^, starting with 
D'^^ = 1, a scalar. To this purpose, we exchange the order of scalars in {<p\t//), 

= ^ M'^' . . . M'^^M'^'* . . . M'^'-*, (101) 
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and transpose again the M-matrices, leading to a hierarchy of bracketed sums, with the sum over 
(Tl innermost. The iteration running from ( - Lio ( - 1 then reads: 



= y M'"D^^^M'"\ (102) 



which can be extended to the calculation of matrix elements as 

^[f-i] = orcKM^iD^(^Nr'\ (103) 

D^"^ then contains the result. 

This approach is very useful for book-keeping, because in DMRG we need operator matrix 
elements for left and right blocks, which is just the content of the C- and D-matrices for blocks A 
and B. As blocks grow iteratively, the above sequence of matrices will be conveniently generated 
along with block growth. 

4.2.2. Transfer operator and correlation structures 

Let us formalize the iterative construction of C''^^ -matrices of the last section a bit more, 
because it is useful for the understanding of the nature of correlations in MPS to introduce a 
transfer (super)operator which is a mapping from operators defined on block A with length 
^ - 1 to operators defined on block A with length i, 

[\a(-i){a'f_^\} ^ [\at){a'f\}, (104) 

and defined as 



'\a-e 

ae-i ,a', , at, a', \ Te 



\ 

{\ae.i){a'^_,\){\at){a'^\), (105) 

(ai.ia'),(ai,a') 



where we read ofl" the expression in brackets as the matrix elements of E^^^ of dimension {D^_^ x 
D^), the M-matrix dimensions at the respective bonds. It generalizes to the contraction with an 
interposed operator at site € as 

= O'^^-'^'MM'"'* ® MM'^^ (106) 

How does E^^^ act? From the explicit notation 

Ef , = yMM-'; .M^'^', (107) 

we can read £'''^[(9''^"^^] as an operation on a matrix d^^'^J'^ as 

£M[0[^-i]] = 2] M'"^0^^-^^M'" (108) 

or on a row vector of length Z)^ j with coefficients Vaa' - O^^'^J^ multiplied from the left, 

38 



The C-matrices of the last section are then related as 



(110) 



but we will now take this result beyond numerical convenience: If Df_i - D(, we can also ask 
for eigenvalues, eigenmatrices and (left or right) eigenvectors interchangeably. In this context 
we obtain the most important property of E, namely that if it is constructed from left-normalized 
matrices A or right-normalized matrices B, all eigenvalues \Ak\ < 1. 

In fact, for Ai - I and left-normalized A-matrices, the associated left eigenvector v^a' - 6aa', 
as can be seen by direct calculation or trivially if we translate it to the identity matrix: 



£[/] = ^A'^' ■/■A'^ = 1 ■/. 



(Ill) 



The right eigenvector for E constructed from left-normalized A matrices is non-trivial, but we 
will ignore it for the moment. For right-normalized B matrices, the situation is reversed: explicit 
calculation shows that Vga' - Saa' is now right eigenvector with A\ - \, and the left eigenvector 
is non-trivial. 



To show that 1 is the largest eigenvalue ||98|], we consider C — E[C]. The idea is that then 
one can show that s'^ < si for the largest singular values of C and C, if E is constructed from 
either left- or right-normalized matrices. This immediately implies that all eigenvalues of E, 
\Ai\ < \. C - AjC implies s[ - \Ai\si, such that \Ai\ > 1 would contradict the previous statement. 
The existence of several = 1 cannot be excluded. The proof runs as follows (here for left- 
normalized matrices): consider the SVD C - U^SV. C is square, hence WU - UW - V^V - 
yy ' = /. We can then write 



C' = ^A'"^f/+5yA'^ = [ (f/Ai)^ ... (UA'')^ ] 



s 




' VA^ 




s 


s 








s 


s 




_ VA'' 




s 



Q. 



(112) 

We have P'P = / and g^g = / (however PP'' + I, QQ'' + /), if the A'^ are left-normalized: 
P'^P - 2^ A""' f/^ f/A"" - A'^^A°" = / and similarly for g; they therefore are reduced basis 
transformations to orthonormal subspaces, hence the largest singular value of C must be less or 
equal to that of S , which is ii . 

Independent of normalization, the overlap calculation becomes 



and expectation value calculations before proper normaUzation by ((AliA) would read 



^IL-2]^[L-1]^[L] 



Ot-i 



(113) 



(114) 



Numerically, this notation naively taken is not very useful, as it implies 0{D^) operations; of 
course, if its internal product structure is accounted for, we return to 0{D^) operations as previ- 
ously discussed. But analytically, it reveals very interesting insights. Let us assume a translation- 
ally invariant state with left-normalized site-independent A-matrices (hence also site-independent 
E) with periodic boundary conditions. Then we obtain in the limit L — > oo 
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Uk 

= Y^{\\E^I^\k)A{-'-\k\E^^\\) (L^cx,) 



where Ak are the eigenvalues of E. We have used \Ak\ < 1 for E from normaHzed matrices 
and that /li = 1 is the only eigenvalue of modulus 1; but relaxing the latter (not necessarily 
true) assumption would only introduce a minor modification. \k) and {k\ are the right and left 
eigenvectors of (non-hermitian) E for eigenvalues A^. (1| corresponds to the eigenoperator / for 
E from left-normalized A. 

The decisive observation is that correlators can be long-ranged (if the matrix elements ( 1 \E^^ \ 1 ) 
are finite) or are a superposition of exponentials with decay length - -1 / In Ak, such that MPS 
two-point correlators take the generic form 



+ Y^Ckt-'l^\ (115) 



k=2 



where r = |; - / - 1| and q = {\\E^^\k){k\E^^\\) for / < 

A simple example can be extracted from the AKLT state of the previous section. The eigen- 
values of E were already found to be 1, -1/3, -1/3, -1/3. For the spin-operators, the matrix 
elements for long-range order vanish, such that the correlation {S'.S'p - (12/9)(-l)-'"'e~^-'"'^'''^ 
for i > i, a purely exponential decay with correlation length ^ = 1/ In 3 = 0.9102. On the other 
hand, for the string correlator, the long-range matrix elements are finite, and long-range order 
emerges in the string correlator 

The form of correlators of MPS has important consequences: the usual form correlators take 
in one dimensional quantum systems in the thermodynamic limit is either the Ornstein-Zernike 
form 

<OoO,)~^ (116) 
or the critical power-law form (maybe with logarithmic corrections), 

<(9oO,)~^-^ (117) 



The AKLT state belongs to a very special state class whose correlation functions mimic a quan- 
tum system in a lower spatial dimension (so-called dimensional reduction), which removes the 
V^-term; the AKLT state sits on a so-called disorder line, where such phenomena occur |99]. 

Any finite-dimensional MPS therefore will only approximate the true correlator by a super- 
position of exponentials. It turns out that this works very well on short distances, even for power 
laws. What one observes numerically is that the true correlation function will be represented ac- 
curately on increasingly long length scales as D is increased. Eventually, the slowest exponential 
decay will survive, turning the correlation into a pure exponential decay with ^ = -I /In A, where 
A is the largest eigenvalue of E that contributes to the correlator. The comparison of curves for 
various D is therefore an excellent tool to gauge the convergence of correlation functions and the 
length scale on which it has been achieved. 
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4.2.3. MPS and reduced density operators 

As we have already seen for DMRG, the concept of reduced density operators is of impor- 
tance in various ways. Let us express them using the MPS notation. We have 

= ^ A"^' . . .A'^^A'^'i* . . .A'^'^*\(t){o-'\ = Yj^"' ■■ • • .A'^i"' |o-)<tr'|. (118) 

(r.CT" (r,(T' 

We now bipartition into AB, where A contains sites 1 through € and B sites { + \ through L. 
Tracing out the degrees of freedom of B in the last expression we obtain 

pf=TrBmil,\= Yj A-K..A''^pfA^^\..A'^^\a-){cr'\, (119) 

(T,(T'€A 

where 

pM = Y A'^'*' . . . A'^'-A'^'-i' . . . A'^'^^K (120) 

This equation inmiediately implies a recursion relation between different reduced density matri- 
ces, namely 

ere 

In the thermodynamic limit L— >oo,/'^ooofa translationally invariant system, we may 
therefore ask whether a fixed point relationship 

pi=j;AVlA-t (122) 

is fulfilled. 

All these equations hold even if the matrices of the MPS are not left-normaUzed. In the case 
that they are, we can directly express the density operator in the orthonormal basis generated by 
the A-matrices, namely 

pf = Y^f)ae4a()AA{a'(\. (123) 

Similar relationships hold for the reduced density operator of B, where (using fi-matrices 
now) we obtain 



pf = Tta\0W = Yj ^'^'■^ ■ ■ ■ B°'''^'^'pfB'^'*' . . . B°'''\(r}{(r'\, (124) 
cr,(r'€B 

where 



(TeA 

and the recursion relationship 

pf = Y B'^'Vb""^'"- (126) 
giving rise to a potential fixed point relationship 

p^ = ;^Z?-tp/B- (127) 

cr 
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Again, all these relationships would hold for arbitrary MPS matrices, but if they are right- 
normalized, we again get an expression in an orthonormal basis, now generated by the B- 
matrices, 

Pb^ = Yj^^B^oe.a^aeyBBia'.l (128) 

In the case of a mixed-canonical state = 2o-^'^' • . ■A'^'^VB"''*' . ..B"''-\(r) a rerun of the 
calculation shows that 

pf = T»P^ (129) 



and 



expressed in an orthonormal basis. 



pf = W^'W, (130) 



4.3. Adding two matrix product states 

An operation that one needs comparatively rarely in practice is the addition of two MPS. Let 
us first consider the PBC case, which is easier. Taking two MPS, with no normalization assumed. 



|(/^) = Yj Tr(M'^' . . . M''^)\a) |^) = ^ Tr(M'^' . . . M'^Olo") (131) 
+ \(P) = Y, Tr(A^' • • • ^^)lo-> (132) 



a (T 

we can write down 



where 

N"^' = M"' ® M"^' . (133) 

This means that we simply take M and M as diagonal blocks of a matrix A^. The diagonal 
block structure imphes that upon multiplying the A' matrices the result is again diagonal, with 
MMMMM ... in the first and MMMMM ... in the second block. Then the trace can be split, 
and we are back at the original states: 

TriNNNNNN) - Tr| ^^^^^^ MmUMM ) " ^'^iMMMMMM)+Tr{M MMMMM). 

(134) 

In the case of OBC, we can proceed in exactly the same fashion. On the first and last sites, 
something special has to happen: naively, the first and last dimensions would go up to 2, and the 
scalar nature be lost. Physically, these indices are dummies anyways. So what we have to do 
(and a simple calculation shows that this works) is to form a row vector [M M] and a column 
vector [M MY on the last sites, from the row and column vectors of the original states. 

Addition of MPS therefore leads to new matrices with dimension Dai = Dm + Dgj, such that 
MPS of a certain dimension are not closed under addition. It is also obvious that in many cases 
this way of describing a new state is uneconomical: the extreme case would be adding |^) + 
where the resulting state is the same, just with a prefactor 2, so matrix dimensions should not 
increase. So after additions it is worthwhile to consider compressing the MPS again to some 
lower dimension, which depending on the states added may or may not (like in the example) 
incur a loss of information. 
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Figure 25: For canonization, sets of matrices on a given site are brought together in a single matrix. 



4.4. Bringing a matrix product state into canonical form 

For a general matrix product state, no particular demands are placed on the matrices M'^' 
except that their dimensions must match appropriately. Certain classes of matrices are to be 
preferred, namely left- and right-normalized matrices, leading to left- and right-canonical MPS: 
certain contractions become trivial, orthonormal reduced bases are generated automatically. 

In order to bring an arbitrary MPS to canonical form we exploit that SVD generates either 
unitary matrices or matrices with orthonormal rows and columns which can be shown to obey 
the left- or right normalization condition. 

4.4.1. Generation of a left-canonical MPS 

Setting out from a general MPS, without normalization assumption, making the contractions 
explicit, 

l*^) = Z Z MZaMZa, ■ • » (135) 



we reshape M^^^ by grouping physical and left (row) index to carry out an SVD on the new M, 
yielding M = A5 : 

O" fl,,... 
(T a],... S] 



(T 02,... S\ 



\ 



\ "1 



= (136) 

cr fl2,... si 

As A'A = / due to SVD, after reshaping to A""', left-normaUzation holds for A°"'. The re- 
maining two matrices of the SVD are multiplied into M°"-, such that a new MPS with M^-^^ - 

Y.atS s^.s^vluat Kla. IS generated. 

Now the procedure can be iterated: M^-^^ is reshaped to M^a-j.so.aj (Fig- US, singular value 
decomposed as AS , generating A(o-, ^j, reshaped to a left-normalized A'^j- ,,. The right two 
matrices of the SVD are again multiplied into the next ansatz matrix, and so forth. After the last 
step, left-normalized matrices A^^j live on all sites. 5'ij(y')i,i, a scalar as A"'^ is a column 
vector, survive at the last site, but this scalar is nothing but the norm of li/*). We may keep it 
separately if we want to work with non-normalized states. 
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This procedure extends trivially to normalization: we identify the prefactor of the state, but 
instead of storing it, we simply set it to 1. As we do not use the singular values explicitly, 
the above procedure can be easily reformulated using QR decompositions, along the lines of 
Section 14.1.31 Standard QR, however, does not show us whether the matrices used are bigger 
than necessary, i.e. have singular values that are zero, such that matrices can be trimmed to a 
smaller size without loss of accuracy; this would only be possible using rank revealing QR; for 
many MPS it is however clear from the underlying physics that the spectrum of singular values 
has a long tail, such that this issue does not arise. The same argumentation holds also for the 
generation of a right-canonical MPS, which we turn to in the following. 

4.4.2. Generation of a right- canonical MPS 

The same procedure can be applied to arrive at a state of right-normalized matrices, by car- 
rying out a sequence of SVDs starting from the right on reshaped matrices {M°") — > M = USB, 
splitting B into matrices Z?°" that are right-normalized (due to BE^ - I), and multiplying U and S 
to the left, creating the matrix to be singular value decomposed next: 



y y ...M°""„ M°""„ M"'' Act) 

/ 1 / 1 01-3,01-1 aL-2,aL-l fli-i,!' ' 

Cr ....at-i 

~ Z Z Z ■ ■ ■ ^01-2,01-2^01-2,01-1 

SL-\,SL-\BsL-l,(crL,l)\0') 

cr ...ClL-\ Sl-\ 

f \ 

= Z Z Z ■ ■ ■ Z 

st-x,SL-\ ^Si-iXo-,,,!)!^) 
O" ...Ol-2 \ol-\ / 

= Y,Y,Y,---Ktio,_Mt-U-K-j^^' (137) 



cr ...az,-2 Sl- 

proceeding as before, with the sole differences that (i) the direction is reversed and (ii) reshaping 
now groups the physical index with the column instead of the row index: M'^l_,,Si ^ ^fli_i,(o-,j,) 
Bsi_i,{o-iSi) — » B'^!-,,sr 

4.5. Approximate compression of an MPS 

The (rare) addition of MPS and various algorithms that can be formulated with MPS lead 
to a substantial increase in matrix dimensions of the result. It is therefore a recurrent issue how 
to approximate a given MPS with matrix dimensions (DJ x D'.^^) by another MPS with matrix 
dimensions (D, x Z),+i), where D, < D'., as closely as possible. 

Fundamentally, two procedures are available, SVD compression and variational compression. 
Both have advantages and disadvantages: for small degrees of compression, D ~ D', SVD is fast, 
but it is never optimal; it becomes very slow if D' » D. Variational compression is optimal, but 
slow if the starting point is chosen randomly, can however be greatly speeded up by providing a 
good trial state e.g. from the SVD approach. Generally, issues of getting stuck in a non-optimal 
compression may arise in the variational ansatz. 

Depending on the specific nature of the state to be compressed, procedures can be optimized, 
for example if the MPS to be compressed is a sum of MPS or if it is the result of the application 
of a matrix product operator (MPO; Sec.|5]l to an MPS. 
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4.5.1. Compressing a matrix product state by SVD 

Let us consider an MPS in mixed canonical representation, 

\,/,) = YjA"''A'^^ . . . A'^'A^fi'^'*' . . . B'^'-'B'"-\o-), (138) 
a 

from which we read off the Schmidt decomposition |^) = Sai\af}A\af}B, where the states on 
A and B form orthonormal sets respectively; this follows from the canonical construction. Let 
us suppose there are D' states each for this decomposition. We now look for the state 10-) that 
approximates ItJ/) best in the 2-norm and can be spanned by D states each in A and B. We have 
shown that SVD provides the result by retaining the D largest singular values, and the compressed 
state simply reads \if/) - Saf\a()jn\ai)B, providing a simple truncation prescription: retain the 
first D columns of A"'', the first D rows of B°''+\ and the first D rows and columns of A'^^'. If 
normaUzation is desired, the retained singular values must be rescaled. 

This procedure rests on the orthonormality of the states on A and B, therefore can only be 
carried out at one bond. In order to shrink all matrices, we have to work our way through all 
mixed canonical representations, say from right to left, truncate, and shift the boundary between 
left- and right-normalized matrices by one site to the left, using techniques from canonization. 

After the first step of right-canonization of a left-canonical state, it reads: 

\^(L-i)y ^Y^^a-, _ __A'^'-'USB^'-\cr), (139) 
cr 

where 1 have already reshaped B, which is right-normalized and guarantees that states formed as 
= Zo-t(S'^'^)aL-i,i|o'L) are orthonormal. But so are the states 

I«l-i)a= Yj iA'^' ■■■A'''-'U\a,.Mi...cTL-i), (140) 

a-i...a-t_i 

as SVD guarantees U^U = 1: we are just doing a basis transformation within the orthonormal 
basis set constructed from the left-normalized A"'' . Hence, we have a correct Schmidt decompo- 
sition as 

= Sa,JaL-l)A\aL-l)B- (141) 

OL-l 

The difference to a right canonization is now the truncation: matrices U, S and B'^'^ are truncated 
(and singular values possibly renormalized) to U, S and B"''^ just as explained before: retain the 
D largest singular values. B"^^ is still right-normalized. The next A"'^-' to the left, U and S are 
multiphed together to form M'^'^-' . By reshaping, SVD and reshaping as 

Kj = Mi,(^j^ = ^ UikSkkBkx^j) = Yj UikSkkKj (142) 

k k 

we obtain right-normahzed B°'^-' , truncate U, S and B"''^-^ to U, S and B"''--^ , and the procedure 
continues: 

10') = ^ A'^' . . . A'^^-M'^^-^ (A'^^-' US ) B"^^ lo-) 
cr 

^ YjA"' ... A'^" A'^" (a'^'^-' us ) B^'-lar) 

(T 
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A°'''-'A'"'-'-M'"'-'B'"-\(t) 



cr 





At the end, the compressed MPS li^) is right-normalized and given by B-matrices. As we will see, 
this compression procedure is just the truncation that is carried out by (time-dependent) DMRG 
or TEBD as they sweep through the chain. Both methods at each bond have correctly normalized 
matrices (i.e. orthonormal states) to the left and right, carry out the cut and proceed. 

The disadvantage of the procedure is that a one-sided interdependence of truncations oc- 
curs: the matrix M always contains a truncated U from the previous step, hence is truncation- 
dependent. Generally, truncations cannot be independent of each other because for each decom- 
position (I138l l truncations of the A- and B-matrices affect the orthonormal systems, but here the 
dependence is one-sided and "unbalanced": truncations further to the left depend on those to 
the right, but not vice versa. If the truncation is small - which it usually is for small time steps 
in time-dependent DMRG - the introduced additional inaccuracy is minor; however, for cases 
where large truncations may occur, the dependence might become too strong and the truncation 
far from optimal. 

A second concern regards efficiency: for matrix dimensions (m x n), m > n, the cost of SVD 
is 0(mn2). This means that the SVD cost 0((D'fdD) if D' < dD and O(D'd^D-) otherwise; 
the matrix multiplications cost 0{dD{D')^). In many applications, D' » D; then this method 
becomes quite slow. The situation is even worse if the original state is not in canonical form and 
has to be brought to that form first by a sequence of SVDs, that are of order (9((D')^). 

Let me conclude this section with the remark that of course we can of course also compress 
by imposing some e which at each truncation we accept as maximal 2-norm distance between 
the original and the compressed state (given by the sum of the squares of the discarded singular 
values), implicitly defining D. 

4.5.2. Compressing a matrix product state iteratively 

The optimal approach is to start from an ansatz MPS of the desired reduced dimension, and to 
minimize its distance to the MPS to be approximated iteratively, i.e. by changing its M"' matrices 
iteratively. The matrices play the role of variational parameters. 

The mathematically precise form of optimal compression of from dimension D' to 
with dimension D is to minimize - \4i)\^, which means that we want to minimize - 
- {4'\4') + (lAliA) with respect to li/^). Let us call the matrices M and M respectively, to 
emphasize that we make no assumption about the normalization. Expressed in the underlying 
matrices M, this is a highly nonlinear optimization problem. 

But this can be done iteratively as follows. Start with an initial guess for which could 
be an SVD-compression of \^), arguably not optimal, but a good starting point. Then we sweep 
through the set of M"'' site by site, keeping all other matrices fixed and choosing the new M"'' , 
such that distance is minimized. The (usually justified hope) is that repeating this sweep through 
the matrices several times will lead to a converged optimal approximation. 

The new M°"' is found by extremizing with respect to Ma'*^ ai^ which only shows up in 
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Figure 26: Linear equation system to be solved for iterative compression of an MPS. The fatter lines correspond to the 
state to be compressed, the thinner lines to the compressed state. The unknown matrix is circled. 



-{ifr\t//) + We find 
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y(M"'*...M'^-'*)i,«,_,(M' 



(T 



. M"^ = 0. 



The sum over cr* runs over all physical sites except /. This system looks compUcated, but is in 
fact quite easy. Keeping the matrix to be found, M°"', explicit, we may rewrite this equation as 



(143) 



If, for each cr,-, we interpret the matrix M°'' as a vector v of length D^, O as a matrix P of 
dimension {D^ x D^) and 0°^' as a vector b of length D^, we have a linear equation system 



Pv = b. 



(144) 



The result v can then taken to be the matrix we are looking for As this system is usually too big 
for a direct solution, an iterative solver has to be used, such as a conjugate gradient method. The 
system is Hermitian, as can be seen from the construction of P: unconjugated and conjugated M 
simply reverse their role under transposition. Once again, the graphical representation is simplest 
(Fig.|26j. 

As in the later case of finding ground states variationally, it is important to realize that the cost 
of the matrix-vector multiplications in the conjugate gradient method is not 0(D*) as dimensions 
would naively suggest. There is an obvious factorization, which becomes particularly obvious 
graphically, that then leads to a cost of O(D^): 



where 



Vo-,-1 



V o-i 



(145) 



(146) 



and similarly R. In the graphical representation (Fig.l28ll. they are simply the contracted objects 
to the left and right of the circled M-matrix we are solving for 
Then 



2 o«,_,«„«;_,«;M^;,,; = X^''--;-. 

47 



(147) 



n 



Figure 27: Equation for iterative compression of an MPS for a suitably normalized state. The fatter lines correspond to 
the state to be compressed, the thinner lines to the compressed state. 



L = 



n 



R = 



Figure 28: Iteratively constructed objects L and R for compression. 



two operations of cost 0{D^). 

A similar decomposition simplifies the calculation of the vector b, which is formed from 



Oalta, as 



2 L,,_,„lM'^l^,.Ra:,a',, (148) 



with L and R as indicated in Fig. [27] In fact, calculating L and R is nothing but carrying out the 
first steps of an overlap calculation, starting from left or right. The result would then be the C 
matrix produced there at intermediate steps. If one sweeps through the system from left to right 
and back one can build L and R iteratively from previous steps, which is the most eflicient way. 

We can however drastically simplify the compression procedure if we exploit the canonical 
form! Assume that \iff} is in mixed canonical form A°"' . . . A""' 'M""'B°"'+' . . . B°'^ and we want to 
update M°"' : to the left of the matrix to be updated, everything is left-normalized, to the right 
everything is right-normalized and the form of M°"' does not matter, as it will be recalculated 
anyways. 

Then L„, j,,,'^ - 6a,_i,a'. ^ because of left normalization, and similarly Rai,a'. - Sai,a'. because of 
right normalization, hence Oa,_ia,,a'. = <5a,_,,fl' Sa,,a'.- In the linear equation system this means 
that P - I and we have trivially 

v^b, (149) 

so there is no need to solve a large equation system (Fig.lZTli. 

To make this work for an entire chain, we have to shift the boundary between the left and right 
normalized matrices as we move through the chain. Assume we begin with all left-normalized 
matrices. Then we move through the chain from right to left, start by solving on the last site for 
M"''-, right-normalize it via SVD (or, as we do not need the singular values, more cheaply by QR) 
as before, to obtain A°"' . . . A""'^ -M""'^ 'B""'^ where M""'^ ' is in general without normaUzation. It is 
now optimized as 



(150) 
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Figure 29: Equation for iterative compression of an MPS in a two-site approacli. 



where 











/ ) 


J 



(151) 



and similarly R, the result is right-normalized, and so on as we go through the chain. At the end, 
all matrices are right-normalized, and we restart from the left. 

In order to assess convergence, we can monitor at each step - and observe the 

convergence of this value; if necessary, D has to be increased. The calculation may seem costly, 
but isn't. If we keep in proper mixed normalization, and use Eq. ( llSOI l to simplify the overlap 
we find 

111^) - ^i-Yj ^m'"'M''% (152) 



which is easy to calculate. The subtracted sum is just at the end, this allows us to normalize 

the state by simple rescaling. 

As already hinted at for single-site DMRG - and we will discuss this issue at length in Sec. |6] 
- there is a danger that this variational ansatz gets stuck in a non-global minimum for the distance 
between the compressed and the original state. Often (but not always) it is helpful to consider two 
sites at the same time, by analogy to two-site DMRG, for optimization. An operation count shows 
that this is somewhat slower Assume the compressed \^) is in a mixed-canonical representation 



(153) 



Running through the same arguments as before, optimizing with respect to M^I_^X+\ ' yields an 
equation as in Fig.|29]for Mat'^^'al+i ■ The major change occurs now: we reshape the new M"''°''+' 
as M(a, ,o-f),(o-,+,flf+,), caiTy out an SVD to obtain 



flfXff+iaf+l) 



(154) 



where the M is formed from reshaping US . In fact, it is discarded because we shift one site 
towards the left, looking for M°''-'°'' . We can also use a cheaper QR decomposition of M' 
instead, to obtain B from g ' ■ 

Let me conclude this section on compression by discussing the relative merits of the methods. 
If the compression is only small, the interdependency of the SVD approach will not matter too 
much. Still, the variational ansatz is superior; its only weakness is that because of its iterative 
nature one has to provide an initial guess for the compressed state. Taken randomly, the method 
will waste a lot of time on just getting into the right vicinity. Therefore, the smart proposal is to 
take the S VD-compressed state as the first input into the iterative method. How can we avoid the 
potentially high costs due to D' at least partially? 

In practice, compressions occur mainly in two situations: (i) MPS have been added, hence 
the matrix dimensions have been added; (ii) a matrix product operator (MPO) has been applied 
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to an MPS; we will see that this leads to the multiplication of the matrix dimensions of MPS and 
MPO. 

In the first case, the variational compression can be speeded up by using the fact that li/r) - 
101 ) + 102) + ■ ■ ■ Then we may rewrite the variational equation as 

If we work out the equation (assuming mixed canonical representation), the right-hand side con- 
sists now of a sum of n overlaps involving D-dimensional matrices, instead of one overlap in- 
volving D- and nD-dimensional matrices (costing up to OirP'D^) and 0{nD^) in the two types 
of matrix multiplications occurring in the overlap) we now have n overlaps costing O(D^). For 
large n, the "decomposed" approach should be up to n times faster. 

The second case we postpone for details until we have discussed MPOs. The idea is to 
carry out an SVD compression, but without the particularly costly step of previously ensuring 
correct normalization; if for some reason the block states are almost orthonormal nevertheless, 
the outcome should be quite reasonable (and can be brought into canonical form, which is ch eap 
after compression) or can at least serve as a reasonable input for the variational method jlOOll . 

4.6. Notations and conversions 

So far, we have explored an MPS notation based on one set of matrices per site; special 
normalization properties for these matrices were exploited to arrive at MPS with attractive ad- 
ditional features (like the generation of orthonormal sets of states or the encoding of a Schmidt 
decomposition). If we consider our lattice with sites 1 through L, it would be useful in view of 
the DMRG construction to be able to access easily all L - 1 possible bipartitionings of the system 
AB that can be obtained with a single cut. 

Such a notation has been introduced by Vidal ll47ll and takes the following form: 

o-i,...,o-i 

where we introduce on each site { a set of d matrices F"''^ and on each bond { one diagonal matrix 
A'^J . The matrices are specified by the demand that for arbitrary 1 < i < L we can read off the 
Schmidt decomposition 

|iA> = ^ Sai\a{}A\at}B (157) 

at 

where the Schmidt coefficients are the diagonal elements of A'^^, - Ajf^|a, and the states on A 
and B are given as 

\ae)A = 2 (F'^'A['Jr'^^..Al^-'JF'^')«.|cri,...,cr,), (158) 

o"i,...,o-f 

\cie)B = (r"'^'A"^'^r-'-...Al^-')F-0«Jcr,+i,...,cri>, (159) 

a-f+i....,CTi 

where the states on A and on B are orthonormal respectively, reminding of similar constructions 
from A- and B-matrices. Graphically, the new notation can be represented as in Fig. [30] It is 
obviously a more explicit version of the A-matrix notation with the advantage of keeping explicit 
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reference to the singular values, reduced density matrix eigenvalues and entanglement: Cutting 
bond {, the reduced density operators pA and pB read in eigenbasis representation 



pM=pM^(AM)2, (160) 

more precisely (but irrelevant for real and diagonal h!^^^) - A'^^A'^^' and p^' = A^^" A^^, 
where the eigenstates of p^^ ^^'^ P^b given by {la^)^} and {la^)^} respectively. The von Neu- 
mann entropy of entanglement can be read off directly from A^^ as Sa\b = -Tr(At^l)^ log2(At^)^. 

Before exploring the connections to other notations, let us first show that any quantum state 
can indeed be brought into that form by a procedure in close analogy to the one that decomposed 
1^) into a product of A-matrices (or B-matrices, for that matter). Starting from coefficients Co-i...o-i, 
we reshape to *Po-i,((rT...(rL)> which is SVDed iteratively. We label the singular value matrices A^. 
After the first SVD, we rename A°"' to r°"' . In the subsequent SVDs, as before we form the next 
matrix to be SVDed by multiplying A and into "P, reshaping such that there is always one 
a- and one cr-index for the rows. Using the reshaping of U(a(_ia-e),ae — > '^ae-uae already used, we 
obtain 



~ ^ r^l'*^(«l0"2),(0-3-0't) 

~ ^ ^al^ai,a2 ^02,02^^^ ^02'(<^^-<^l) 



C(a20-3),(0-4...0-i) 



ai ai,a\ 01,02 M 
01,02 

ZT^ia[1] 7^2 40-3 A [3] /yU 
^ oi ^^^01,01 01,02-^02,03 ■"■03,03^' )ai,(a'i-o-L) 



- z 



pO-i.[l] W7-2 »[2] po-3 m 

^ 01 ■''^01,01'^ 01,02 ■"■02,02 02,03 ^ (asa-Mo-s.-CTL) 



and so on. The crucial difference to the decomposition into A-matrices is that each A is decom- 
posed, using the knowledge of A^^"^^ obtained in the previous step, into 

K' a =^a~^i o. (161) 

Of_i,0^ 0^_l,Of_i Of_i,Of' ^ / 

which implies a division by the diagonal elements of A^^"''. If in our SVD we keep only the 
non-zero singular values, this is a mathematically valid operation, albeit potentially fraught with 
numerical difiiculty. Ignoring this issue in this conceptual demonstration, we do arrive at a 
decomposition of the desired form; in order to prove that it is indeed correct, we have to show that 
at each iteration we indeed obtain a Schmidt decomposition. But this is easy to see: The matrices 
to the left of any A^^^ can aU be grouped into (or rather, have been generated from) left-normaUzed 
A-matrices, which generate a set of orthonormal states on the part of the lattice ranging from site 
1 to {. On the right hand side of any A^^^, there is a matrix with orthonormal rows, which 
means that the states la^)^ = Xa-e+,,...iV^)a{,a-ui - \'^e+i • • •) are also orthonormal. Hence, the SVD 
giving A^^ indeed leads to a valid Schmidt decomposition. 
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r°i Afi] r°2 A[2] r"L-iA[L-u 




A A A A B B B B B 

Figure 30: Representation of an MPS in Vidal's notation. Singular values remain explicit on bonds (diamonds). A sit 
on bonds, T on sites. By construction, adjacent A and T can be contracted to A or B matrices, that are either left- or 
right-normalized. The state can be trivially grouped into a string of A (giving orthonormal block states), a singular value 
matrix, and a string of B (giving orthonormal block states). 



An alternative way of obtaining this notation would be to carry out a standard left-normalized 
decomposition, and store all singular value matrices generated (and previously discarded) as A''^, 
and to insert afterwards the identities A['^(A''1)"' between all neighbouring A-matrices A""' and 
A°"'+' . Then using Eq. (II6II 1 leads to the same result. 

Similarly, starting the decomposition from the right using the right-normalization of B- 
matrices the same state is obtained with a grouping 

Bn' „ = A[^, , (162) 

where for notational simplification for this and for the corresponding equation for the A-matrix, 
Eq. (1161b . it is useful to introduce dummies A'"' and A'^^ that are both scalar 1 . 

The groupings for A and B-matrices allow to reexpress the left- and right-normaUzation con- 
ditions in the FA-language: The left-normalization condition reads 

I^Yj = 2 r^^' ' Al'-'l^Al'-'lF'"' (163) 

or, more compactly, 

2r-'Vl;"'Jr-' =/. (164) 

The right-normalization condition reads 

^r>wr-,t^;. (165) 

cr, 

Interestingly, Eqns. ( 1164b and ( 1165b also arise if we translate the density operator recursions 
Eqns. (1121b and ( 1126b using Eqns. ( 1161b and ( 1162b . A matrix product state in the form of 
Eq. (1156b which meets the constraints Eq. ( 1164b and Eq. (1165b is called canonical. 

Conversions between the AB-notation, the FA-notation and also the block-site notation of 
DMRG are possible, albeit fraught with some numerical pitfalls. 

Conversion FA — » A,B: The conversion from FA — » A,B is easy. If one introduces an 
additional dummy scalar A''^^ = 1 as a "matrix" to the very left of we can use the above 
defining Eq. ( 1161b to group 

(A™F'^')(A"^F'^'-)(A[2¥'^') . . . ^ A'"A°''-A'^' . . . (166) 
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r°i A[i] r°2 A[2] 




Figure 31: Vidal's MPS notation. A, B-matrix MPS notation, and DMRG block notation. The A-matrices generate the 
left block states, the B matrices generate the right block states. The matrix A''^' connects them via singular values. 
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Figure 32: Representation of a state in single-site DMRG: translating Vidal's MPS notation and A, B-matrix MPS notation 
into DMRG block notation. The A-matrices generate the left block states, the B matrices generate the right block states. 
The matrix elements of 'P°"f are just the coefficients of the DMRG state. 



or using Eq. ( I162l l 



(167) 



In view of what DMRG and other MPS methods actually do, it is interesting to consider 
mixed conversions. Consider bond {! between sites { and i + \. We could contract AT — > A 
starting from the left, giving left-normalized matrices, and FA — > B from the right, giving right 
normalized matrices, leaving out just A''^^ in the center (Fig.lSTTi: 



(A[0]r)(A[ilr) . . . (A[^-2]p)(^K-i]p)^[f](p^K+i])(p^[f+2]) (p^[L])^ 



(168) 



As the bracketing to the left of bond ( generates left-normalized A-matrices and right-normalized 
matrices B on the right, we can multiply them out as done in the recursions of the previous 
Section to arrive at orthonormal block bases for A and B, hence at a Schmidt decomposition 



(169) 
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r°i AH] AP] r°'i-iA[^-i] t"^ 




A A A At^-i] r AM r A[^'+i] B B B 




\a{i-\)A \'^f) \0(m) \ae.+\)B 

Figure 33: Representation of a state in two-site DMRG: translating Vidal's MPS notation and A, B-matrix MPS notation 
into DMRG block notation. The A-matrices generate the left block states, the B matrices generate the right block states. 
The elements of matrix m'^t'^M are just the coefficients of the DMRG state. 



What is more, we can also take one site {€) or two sites {(.,{+ 1) and multiply all matrices into 
one there (Fig.|32]and Fig.l33Ti: 

(A™r)(At'ir) . . . (At^"2ir)(A[^"'irA[^])(rA[^+'J)(rA[^+2J) . . . (rA'^i). (i70) 

Calling the central matrix ^"'i - A^'^^'^F'^'^A''^^, we can write 

\t/f} = Y^A'^' . . .A'^'-'^V'^'B'^'*' . . .B'^'lo-) (171) 
cr 

or, again building block bases, 

l«A>= ^''c-i^AXi.a.laeyB- (172) 

If we group even two sites, we have 

(A™r)(Al'lF) . . . (A[^-2]p^(^[f^i]p^[f]p^[^+i])(p^[f+2]) (p^M) 

or, with central matrix = At^-ilp'^f Alflr'^r.i aI^+U^ 

= ^ A"^' . . .A'^f-'»P'^f'^f+'B'^" . . .fi'^^lo") (174) 
cr 

or, using block bases, 

af-l,flf+l,o"f,o-f+i 

These are just the states considered by "single-site" and the original "two-site" DMRG, which 
keep one or two sites explicit between two blocks. 

Conversion A,B ^ FA: Conversion in the other direction A, Z? — > FA is more involved. 
The idea is to obtain iteratively the Schmidt decompositions (hence singular values) of a state, 
hence the diagonal matrices A''^^, from which the r°"-matrices can be calculated. The procedure 
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Figure 34: Conversions of representations: c is tlie explicit representation by the exponentially large number of state 
coefficients; A, B and FA stand for left-canonical, right-canonical and canonical MPS; M stands for an arbitrary MPS. 
Solid lines indicate computationally feasible conversions, dashed lines more hypothetical ones. 

is in fact very similar to that used to show the existence of a canonical FA representation for an 
arbitrary quantum state. 

Let us assume that is right-normalized, then state coefficients take the form Z?°"' B"'- B"'' B"'" . . 
Then we can proceed by a sequence of S VDs as 

= VM'^'-B'^'B'^* ... 

= r°''{A'^'-A^^^V"')B'^'B°'' ... 

- r°"i A''^r°"w°"-'B°"'' . . . 

and so forth. Here, the to-be-SVDed matrices M""' - A^^^^W^ B°'' . The F""' -matrices are obtained 
from the A""^ -matrices by remembering At'^"'' and using Eq. ( 11611 ), which implies a division by 
singular values. 

The division by singular values is a numerical headache as they can and will often be very 
small, in particular if a high-precision calculation is attempted and even very small singular 
values will be carried along. It is numerically advisable to proceed as in the calculation of the 
(pseudo)inverse of an almost singular matrix and set all < e, with, say, e = 10"^, to and 
exclude them from all sums (e.g. in a Schmidt decomposition). As we order Sa by size, this 
implies shrinking matrices U and V ' accordingly. These small singular values carry little weight 
in the reduced density operators (their square), hence the loss of accuracy in the state description 
is very small compared to the numerical pitfalls. In fact, the problem is that at various places in 
algorithms we implicitly rely on (ortho)normality assumptions that may no longer hold after a 
"wild" division. 

Let me conclude this long section by summarizing the various conversion and canonization 
procedures in a diagram (Fig. [34b. where it should however be kept in mind that some conversions 
are only possible in theory, not in numerical practice. 
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5. Matrix product operators (MPO) 



If we consider a single coefficient {cr\\li) of an MPS, 
it is a natural generalization to try to write coefficients {cr\0\cr') of operators as 151 



{(T\d\(r') = W'"^ W-"^- . . . W'-"''- (176) 

where the W"''^ are matrices just like the M°", with the only difference that as representations of 
operators they need both outgoing and ingoing physical states: 



6-Z 



^<r,<r, ^<r,,T, ^ ^ ^ ^tr,_„T,_, jy'^i'^L |o-)(o-' |, (177) 

cr,(T' 



with the same extension to periodic boundary conditions as for MPS. The pictorial representation 
introduced for MPS can be extended in a straightforward fashion: instead of one vertical line for 
the physical state in the representation of M, we now have two vertical lines, one down, one up, 
for the ingoing and outgoing physical state in W (Fig. [35]). The complete MPO itself then looks as 
in Figure [36] If we want to use good quantum numbers, the methods for MPS translate directly: 
we introduce an ingoing local state quantum number from the top, an outgoing one towards the 
bottom, and an ingoing quantum number from the left and an outgoing one to the right. The rule 
is, as for MPS, that the total sum of ingoing and outgoing quantum numbers must be equal, or 
M{\ai)) + M{\bi-.i)) - M{\o-'.)) + M{\bi)), where I have interpreted the bond labels as states for the 
notation. We may also think about dummy indices before the first and after the last site as in an 
MPS, which reflect in which (definite!) way the operator changes the total quantum number. For 
a Hamiltonian, which commutes with the corresponding operator, the change is zero, and we can 
ignore the dummies. The MPOs we are going to build can all be shown to have good quantum 
numbers on the bonds, because they originate either from SVDs (e.g. for time evolutions) or 
from rules that involve operators with well-defined changes of quantum numbers (e.g. for MPOs 
for Hamiltonians). 

In fact, any operator can be brought into the form of Eq. ( 11771 ). because it can be written as 

O ^ ^ C(a-i...a-L)(a-\...a-'i^)\0-U...,Cri){cr\,...,0-'^\ 

cr\,...,a-L,a-\,...,a-'^ 

C(^io-;)...(o-iO-y|o"i,...,cri)(cr'j,...,cr^| (178) 

cr,,...,o-i,o-'|,...,o-^ 

and we can decompose it like we did for an MPS, with the double index a-jcr'. taking the role of 
the index cr, in an MPS. 

As for MPS, we have to ask how we operate with them and how they can be constructed in 
practice, because the naive decomposition might be exponentially complex. As it turns out, most 
operations run in perfect analogy to the MPS case. 

5.7. Applying an MPO to an MPS 

The application of a matrix product operator to a matrix product state runs as 



cr,cr' 
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Oi Ol 

^ b\ b(_i ^ bf bL.i ^ 

a'l o'ji o'l 

(i) (ii) (iii) 

Figure 35: Elements of a matrix product operator: (i) a corner matrix operator wj ' 'at the left end of the chain; (ii) a 

bulk matrix operator wjj^' ^ 'i^' ; (iii) a comer operator wj^^^'l^'j'^'' at the right end: the physical indices points up and down, 
the matrix indices are represented by horizontal lines. 

Ol Oi 

H H M t 

Ol Ol 

Figure 36: A matrix product operator acting on an entire chain: the horizontal matrix indices are contracted, and the 
MPO is ready to be applied to an MPS by simple contraction of vertical (physical) indices. 

o- a.b 
o- a,b 

o- a,b 

cr 

The beauty of an MPO is that it leaves the form of the MPS invariant, at the prize of an 
increase in matrix size: the new MPS dimension is the product of that of the original MPS and 
that of the MPO (Fig.O. 

The result can be summarized as \(p} = 0\ifr} with |0) an MPS built from matrices N°'' with 

cr'. 

If we use (additive) good quantum numbers, one can show from the sum rules at each tensor that 
they are additive on the in- and outgoing horizontal bonds. 

Once again, a seemingly exponentially complex operation (sum over exponentially many cr) 
is reduced to a low-cost operation: the operational count is of order LcfD^D~, Dw being the 
dimension of the MPO. 
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Figure 37: A matrix product operator acting on a matrix product state: matcliing pliysical (vertical) indices are contracted, 
a new matrix product state emerges, with multiplied matrix dimensions and product structure in its matrices. 

5.2. Adding and multiplying MPOs 

Operations with MPOs follow very much the lines of MPS. If we consider the addition of 
two operators, O and P, that have MPO representations W"'''^' and W'^'°''i , then the resulting MPO 
is formed exactly as in the case of MPS, by the direct sum ^^'^> for all sites 1 < / < Z., 
with the same special rules for sites 1 and L. In essence, we (again) just have to consider cr, and 
cr'. as one "big" physical index. 

The multiplication of (or rather, subsequent operation with) two operators, PO, can also be 
understood easily: the application of O to some state li/r) leads to a new MPS with matrices 
■^(fc ,a i)(fcfl) ~ 2o-' ^'b''^^b-^'ai-\ar Then the subsequent operation of P gives a new MPS with 
K"^-' - - Yjrr' W'Z'°''~ nZ' ,. But from this we can read off right away (and this 

is also obvious from the graphical representation of a sequence of two MPOs applied to a state) 
that the new MPO (with matrices Y°''°''i) is given by 

Hence, MPO dimensions simply multiply as for tensors. If we consider an MPS as an MPO with 
dummy indices in one physical direction, the rule for applying an MPO to an MPS follow as a 
special case. 

5.3. Compressing MPOs and MPO-MPS products 

As for MPS, the question of compressing an MPO may arise. This should be obvious from 
the last section, where MPO dimensions summed up or multiplied. If it is no option to shift the 
issue of compression to the application of an MPO to an MPS (then of dimension DwD and a 
natural candidate for compression), we have to compress the MPO. A typical example would be 
given by the representation of a longer-ranged Hamiltonian in MPO form, which quickly leads 
to large dimensions. 

But we can apply the same techniques as for compressing MPS, both by S VD and iteratively, 
in order to approximate the exact MPO by one with smaller Dw- The only change is that instead 
of one physical index cr, we have now two physical indices cr, cr', which we may take as one 
single index (cr, cr'). The approximation is then done in the Frobenius norm, which naturally 
extends the 2-norm of vectors we used for MPS approximations. 
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At this point it is worthwhile mentioning that it has been proposed that in the special 
case of compressing an MPO-MPS product, an important speedup over the standard methods 
may be achieved: SVD may be very slow if normalization has to be carried out first at a cost 
O(D^D^), but a good starting point for the variational method would be essential to have. But 
the proposed solution from SVD compression may not be bad if the block states are almost 
orthonormal and it seems that in the MPO-MPS product case this is essentially true if both the 
MPO and the MPS were in canonical form (for the MPO again formed by looking at the double 
index as one big index), which can be achieved at much lower cost iOidD^) and OiSo]^), where 
Djv <K D usually, versus O(dD^D^)). Even if the proposed compression is not too good, it will 
still present a much better starting point for the variational compression. So the procedure would 
be: (i) bring both MPO and MPS in the same canonical form; (ii) do SVD compression, of course 
only multiplying out MPO and MPS matrices on the fly; (iii) use this as variational input if you 
don't trust the result too much. 



6. Ground state calculations with MPS 

Assume we want to find the ground state of some Hamiltonian H. In order to find the optimal 
approximation to it, we have to find the MPS \tfr} of some dimension D that minimizes 

E^)ILJ1L_ (181) 

im 

The most efficient way of doing this (in particular compared to an imaginary time evolution 
starting from some random state, which is also possible) is a variational search in the MPS 
space. In order to make this algorithm transparent, let us first express H as an MPO. 

6.1. MPO representation of Hamiltonians 

Due to the product structure inherent in the MPO representation, it might seem a hopeless 
task - despite its guaranteed existence - to explicitly construct a compact MPO representation 
for a Hamiltonian such as 

i-i J J 

^ = Z 2^'^-i + ^^^ti+ J'S]Sl,-hY,S]. (182) 
(=1 (■ 

This common notation is of course an abbreviation for sums of tensor products of operators: 

H - J~S\®S\®I®I®I...+ 



It is however surprisingly easy to express this sum of tensor products in MPO form ['70'] - to 
this purpose it is convenient to reconsider the building block W^'^T combined with its associated 
projector \a-){cr'\ to become an operator-valued matrix Wbb' - Eo-o-' ^ht' l''')^''''! ^^^^ that the 
MPO takes the simple form 



59 



Each W^'^ acts on a different local Hilbert space at site i, whose tensor product gives the global 
Hilbert space. Multiplying such operator-valued matrices yields sums of tensor products of op- 
erators such that expressing in a compact form seems feasible. 

To understand the construction, we move through an arbitrary operator string appearing in H: 
starting from the right end, the string contains unit operators, until at one point we encounter one 
of (in our example) 4 non-trivial operators. For the field operator, the string part further to the left 
may only contain unit operators; for the interaction operators, the complementary operator must 
follow immediately to complete the interaction term, to be continued by unit operators further to 
the left. For book-keeping, we introduce 5 corresponding states of the string at some given bond: 
state 1: only units to the right, states 2,3,4: one 5^, S^, S~ just to the right, state 5: completed 
interaction or field term -hS ' somewhere to the right. Comparing the state of a string left and 
right of one site, only a few transitions are allowed: 1 — > 1 by the unit operator /, 1 ^ 2 by 5^, 
1 ^ 3 by 5", 1 ^ 4 by 5M ^ 5 by -hS\ Furthermore 2 ^ 5 by (7/2)5 3 ^ 5 by (7/2)5 + 
and 4 — > 5 by J~S to complete the interaction term, and 5 — > 5 for a completed interaction by 
the unit operator /. Furthermore all string states must start at 1 to the right of the last site and 
end at 5 (i.e. the dimension Dw of the MPO to be) to the left of the first site. This can now be 
encoded by the following operator- valued matrices: 



and on the first and last sites 
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-hS' 


(7/2)5- 


(7/2)5 + 


7=5- 
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(184) 



W^^^^[-hS~ (7/2)5- (7/2)5+ 7=5= /] 



/ 

5 + 
5- 
5 = 
-hS' 



(185) 



Indeed, a simple multiplication shows how the Hamiltonian emerges. Inserting the explicit op- 
erator representations gives the desired W"''^ -matrices for the MPO. It is therefore possible to 
express Hamiltonians exactly in a very compact MPO form; a similar set of rules leading to the 
same result has been given by | lOlJ- 

For longer-ranged Hamiltonians, further "intermediate states" have to be introduced. Let us 
consider a model with just 5 =5=-interactions, but between nearest and next-nearest neighbours. 



H 



i+2- 



(186) 



Then the bulk operator would read 
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(187) 
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While the 7i -interaction can be encoded as before (moving as 1 ^ 2 ^ 4), for the next-nearest 
neighbour interaction, one has to insert an additional step between 2 and 4, an intermediate state 
3, where exactly one identity is inserted (moving as 1 — > 2 — > 3 — > 4). It merely serves as a book- 
keeping device. Similarly, one can construct longer-ranged interactions. Except the top-left and 
botton-right corner, the non-vanishing parts of W^'^ are all below the diagonal by construction. 

It might seem that for longer-ranged interactions the dimension Dw will grow rapidly as more 
and more intermediate states are needed (one additional state per unit of interaction range and 
per interaction term). While this is true in gene ral, important exceptions are known which can 
be formulated much more compactly 11741 llOlll ; consider for example the following exponen- 
tially decaying interaction strength J(r) - Je^''^^ - JA'^, where r > and A - exp(-l/^). An 
interaction term 2, Ji'')^]^]^r would be captured by a bulk operator 



7 

5-^ ^7 
JAS'~ I 



(188) 



But even if such a simplification does not occur, it turns out that MPOs with quite small 
dimensions and moderate loss of accuracy can be found, either by approxim ating an arbitrary 
interaction function J{r) by a sum of exponentials coded as above 1721 llOlh . minimizing the 
L2 distance \\J{r) — Q'/^' ll in Oi,Ai, where n is given by the Dw and loss of accuracy one 
is willing to consider. Alternatively |74], one can of course construct the exact MPO where 
feasible and compress it by adapting MPS compression techniques to an acceptable Dw (and 
loss of accuracy). 



6.2. Applying a Hamiltonian MPO to a mixed canonical state 

Let us consider in the following mixed canonical representation, identical to the single- 
site DMRG representation. 



or 



(T 



(189) 



(190) 



Let us now look at the matrix elements (a^_icrfflf|77|fl^ jCr^a'^) obtained using the MPO repre- 
sentation for 77. By inserting twice the identity 7 = Y^cr we obtain (the sums with a star 
exclude site () 



{a(-iO-eac\H\a'i_^o-'^a'i) 
= YjYj ■ ■ ■ '^'"'■'"''{ac-icria(\o-){(T'\a'^_y,a'^) 



(T (T' 



<flf_i|cri . . .cr(_i){a[\cri+i . ..o-i}{o-\ . . . cr'^_j|aj,_j)<cr^^, . . .cr^lfl'^) 



(Tt cr'. 
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Figure 38: Representation of tlie DMRG expression {af-i(Tca(\H\a'^ ^(T'^a'^) in MPO/MPS language. Tlie Hamiltonian 
MPO is contracted witli four block state expansions in MPS forni (two bras, two kets, two on block A, two on block B). 
The contracted network decouples into parts L, W and R, con'esponding to blocks A and B and the center site. 



(A^- . . .A'^'-' (B'"'*' . ..B''%JA< . . .A'^^-.)i,«;_,(B"'^-' • ■ .B"^)a^,i 



y y A^'^wri'^'Af', yAr. <t^a"^ 

/ i / I I.fl] l,b[ l,fl, /_i "i'«2 b,,b2 a\,a\ 
\ai,bi,a'.\\a-\a^^ jyo-ijf^ 

\ I 

y B--* wr"""'^'B"h ... y B'^'* , 

flf.af+i b[,b[^\ "f.^f+l i—i "i-i.l «l_pl 



All the beauty of the MPO formulation seems gone, but a graphical representation restores it 
(Fig. [38]). It can be understood most easily from the second or third line of the explicit expressions 
above: the Hamilton MPO (expressed in the product basis) is projected on the block states of A 
and B, which have an expansion in the cr-basis. 

In fact, we can also encode the obvious tripartite structure of the expression as 



(191) 



hi-i,br 

where L and R contain the contracted left and right parts of the graphical network 



Ri 



{aj,bj,a'.;i<{-\\ 



y y A7'-wrr"''Af', ... y a^'--* w^^f^A^r , ;i92) 

/ 1 / 1 l.fli 1,^1 l,a\ /_i at-2,ac-\ b[-2,bl-l <'f-2'''t-l 



y y b:^;'* wr7'"^*'B"h ... yB^^^-.^r^'^fB"^ , 

{ai,bi,a'.;i>{)\crr+iiT'i^^ 



(193) 



We can now write the action of on a state \t//} in the mixed canonical or single-site DMRG 
representation as 

As we will discuss in an instant, H\t//} is the key operation in an iterative ground state search. 
Evaluating this expression naively is inacceptably slow; it can be drastically accelerated on two 
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Figure 39: Update from fl'^'l to F''' by contracting with A'''*, W''' and A'''. While it makes sense mathematically to 
consider the three added tensors as one object, in numerical practice, they are contracted into the network sequentially 
for efficiency. 



counts: first, L and R can be built iteratively in order to maximally reuse available information; 
this involves an optimal arrangement of a network contraction. Moreover, the final action of L, 
R and W on can also be arranged highly efficiently. 

Let us first consider building L and R. In actual applications, we will never carry out the full 
network contraction that stands behind them, because in the spirit of DMRG we are looking at 
blocks that are growing and shrinking in size site by site. The construction of L and R, however, 
is iterative in a way that directly matches block growth and shrinkage. I will illustrate it for L, 
using A-matrices; left-normalization will be exploited explicitly for further simplification at one 
point only such that the formulae are generic. We start by considering the block of size 1 : we 
contract Al^'l and A''^ ' with W''^ . The block basis representation is then given by 



711] 



a-i,a-'.,ao,ho,a'„ 



hoM ^ 



'"0 ao,bo,a' a' a'. 



(195) 



where we have introduced a dummy scalar F^^^j^^ ^, = 1, and where aq, bo, aj, can just take the 
value 1 ; this is just to make the first step more consistent with all that follow. The resulting object 
is a tensor f^^'^^j , corresponding to the three legs sticking out. 

We can now simply continue to contract A, A^ and W on the next site, and the contraction 
update reads 



2 



bi-i,b, 



[i-l] 



and can be represented pictorially as in Fig. [39] 

This construction can be calculated most efficiently by optimal bracketing as 



(196) 



-I] 



, AT' , 

a. , .a. 



(197) 



Here, we have contracted the three new tensors into the network one by one, at operational 
counts O(dD^Dw) in the innermost bracket, then OicfiD^D^) and last O(dD^Dw)- In fact, the 
second operation is faster in practice, as we know that most operators in W are simply zero; the 
remaining ones also often have a simple structure. Another acceleration is possible in the case 
of building L from left-normalized matrices for indices = Dw, if we build H following the 
rules outlined in the previous section: we know that in this case only identities operate towards 
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the left, implying that F^^^[,^^ „> - 5„^,„', simplifying both the innermost bracket and the outermost 
operation. The same idea applies for indices bi - 1 on the right side for building R from right- 
normalized matrices. 

Note that this construction is a generalization of the representation update explained for 
DMRG: a typical situation is the representation of OiOj, where the two operators act locally on 
sites i and j respectively. Then the MPOs are of dimension (1x1) everywhere and W"'"' = da;a-' 
everywhere but on sites / and j, where they read W'^'"'' = O^'^"'"'^' and similarly for j. Pushing 
forward the contractions, F^' is still a scalar 1. Then 

= Ol'l'^'-'^At'l'^'UW'^ (198) 

where F^^\ ^[!]o"it ^ud A^'^"^' are matrices and a multiplication A^'^"''^A^'^'^t is implied. The F^'^- 
matrix is just the operator representation in the block basis, comprising sites 1 through i. 
The update up to site j - 1 then simplifies to 

= 2] A^'^'VI'^-'UW'^*, ii<k< j) (199) 

matrix multiplications implied, and at site j we get again a non-trivial step, 

Ft^T = Yj 0^'^'''''^JA^J'^"''''F^j-^'^A^'^'^i, (200) 

after which updates continue as on the previous sites. Making the matrix multiplications explicit, 
one sees that this is just the construction discussed for the DMRG algorithm. 
In the end, H\>//) can be bracketed advantageously as follows: 

^^ 

\ae-i)A\o-{)\a{)B, (201) 



which scales at worst as 0{D^). More precisely, the innermost operation is 0(D^Dwdy, the next 
one is O(D^D^d^), after this we have a sum of cost 0{D^D^cf). It is advantageous to keep track 
of the structure of W, namely exploiting for which {be-i,bt) configurations it is zero and nothing 
has to be calculated (usually, for most of them), and to use the simplifications for L and R just 
discussed if the state is in mixed-canonical form. 

6.3. Iterative ground state search 

Let us now turn to the algorithm. Assume H given in MPO form and consider a class of 
MPS with predefined matrix dimensions (simply think about a random MPS with matrices M"' 
of desired shape and size, but no normalization assumed for the moment). In order to find the 
optimal approximation to the ground state within this class, we have to find the MPS that 
minimizes 

E=<f^. (202, 

WW) 

It turns out that this can be turned into a ground state algorithm much more efficient than imag- 
inary time evolution from some random state. In order to solve this problem, we introduce a 
Lagrangian multiplier A, and extremize 

- (203) 
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Figure 40: Network to be contracted to obtain the functional to be extremized to find the ground state and its energy. The 
left-hand side represents the term (i//\H\i//), the right-hand side the squared norm 



in the end, will be the desired ground state and A the ground state energy. The MPS network 
that represents Eq. ( I203l l is shown in Fig.l40l 

The problem with this approach is that the variables (the matrix elements M'^^,) appear in 
the form of products, making this a highly non-linear optimization problem. But it can be done 
iteratively, too, and this is the idea that also drives DMRG: while keeping the matrices on all 
sites but one {£) constant, consider only the matrix entries Mgl_^„^ on site £ as variables. Then 
the variables appear in Eq. ( 1203b only in quadratic form, for which the determination of the 
extremum is a benign linear algebra problem. This will lower the energy, and find a variationally 
better state, but of course not the optimal one. Now one continues to vary the matrix elements 
on another site for finding a state again lower in energy, moving through all sites multiple times, 
until the energy does not improve anymore. 

Let us first consider the calculation of the overlap, while keeping the chosen M°'' explicit. 
We find 

(m = Z Z Z K..i,Kl.Kly,Ka', (204) 

where 

K-ual, = Z )«,_,,«;_, (205) 

o"i,...,o-f_i 

K,a', = Z (M'^''' ■■■M'^'M'^'^' ■■■M'"*'%,ar (206) 

o"f+l,--,o"t 

As is particularly clear in the graphical representation, for obtaining the last two expressions the 
same rules about smart contracting apply as for overlaps; moreover, if we move through sites { 
from neighbour to neighbour, they can be updated iteratively, minimizing computational cost. In 
the case where sites 1 through £- 1 are left-normalized and sites f+l through L right-normaUzed, 
normalization conditions lead to a further simplification, namely 

^a[_l,a'^_^ ~ ^af-i,a'f_j ^ a[a[ ~ ^aia'^- (207) 

Let us now consider {tfr\H\i(/}, with H in MPO language. Taking into account the analysis of 
H\iff} in the last section, we can immediately write 

<^i^i'A> = z z z z cr"<:i<'"^<-:..<,.; (^os) 

with L and R as defined before; how such an expression can be evaluated efliciently has been 
discussed previously. 
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Figure 4 1 : Generalized eigenvalue problem for the optimization of M ,a[ ■ The unknown matrix is circled on the left 
and right networks. 

If we now take the extremum of Eq. (I203l l with respect to M^/' we find 

y y y Lr'-^'^-'wr'-iC'^'^M"^ , ,m7 , =o. (209) 

/ 1 / 1 / 1 Or-] b[-i,hc be "r-\'"i ' ' "t-\,a^^ aia^ "i-\'"r 

This is in fact a very simple eigenvalue equation; if we introduce matrices H and by reshaping 

fj _ V r^f-i'^'f-i Ti/°"f'°"f jfi-'''r „„j AT _ vuA \t)B r „„ 

n(o-,a,_,flf),(o-;a;_,cip - L,br-,,br ^bf_i b,_i,b[ bf "i"" ■'^(o-faf_,a,),(o-^f4_,fl'f) - ^ af_ua'f_j at,a'^"°'iK 

well as a vector v with Va-fac^iae = -^afii.flfj we arrive at a generalized eigenvalue problem of 
matrix dimension (dD^ x t/D^), 

Hv-ANv^O, (210) 

represented in Fig. [41] Solving for the lowest eigenvalue Aq gives us a v"^^^^^^, which is reshaped 
back to M^'^a,, Aq being the current ground state energy estimate. 
A few remarks are in order 

• The problem is Hermitian; both H and are Hermitian, as can be seen from the construc- 
tion and the Hermiticity of the MPO employed. 

• In general, dD^ is too large for an exact diagonalization, but as we are only interested in 
the lowest eigenvalue and eigenstate, an iterative eigensolver that aims for the ends of the 
spectrum will do. Typical methods are the Lanczos or Jacobi-Davidson large sparse matrix 
solvers. The speed of convergence of such methods ultimately rests on the quality of the 
initial starting or guess vector As this eigenproblem is part of an iterative approach to the 
ground state, the current M°'' is a valid guess that will dramatically speed up calculations 
close to convergence. 

• Generalised eigenvalue problems can be numerically very demanding, if the condition 
number of becomes bad. But this is no issue for open boundary conditions, if one 
ensures that the state is left-normalized up to site (—1 and right-normalized from site ^ -i- 1 
onwards. Then the simplifications for 'I''^ and imply that is just the identity matrix 
/. The eigenvalue problem then simplifies to a standard one, 

yy y c''''-<'i/?r''M"^ ,-am:' ,^o. (21 d 

/ 1 / 1 / 1 be-] bc-i,b[ be "r-y'tf a,-],ai \ i 

cr\ a\ ^a\b,-],b, 

or Hv - Av - 0, as represented in Fig.|42] The evaluation of the sums will be done using 
the optimal bracketing for To achieve this simplification, one will sweep the position 
t from right to left and vice versa through the chain, such that the optimal normalization 
configuration can be maintained by a single step of the left or right canonization procedure 
after each minimization. 
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Figure 42: Standard eigenvalue problem for the optimization of M ^f ^ ai- The unknown matrix is circled on the left 
network. 

The optimal algorithm then runs as follows. 

• Start from some initial guess for which is right-normalized, i.e. consists of Z?-matrices 
only. 

• Calculate the /^-expressions iteratively for all site positions L - 1 through 1 iteratively. 

• Right sweep: Starting from site { - \ through site L - 1, sweep through the lattice to 
the right as follows: solve the standard eigenproblem by an iterative eigensolver for M°"', 
taking its current value as starting point. Once the solution is obtained, left-normalize M°'' 
into A°'' by SVD (or QR) to maintain the desired normalization structure. The remaining 
matrices of the SVD are multiplied to the M°"'+' to the right, which will be the starting 
guess for the eigensolver for the next site. Build iteratively the L expression by adding one 
more site. Move on by one site, { ^ ( + \, and repeat. 

• Left sweep: Starting from site t - L through site 2, sweep through the lattice to the left 
as follows: solve the standard eigenproblem by an iterative eigensolver for M"'' , taking its 
current value as starting point. Once the solution is obtained, right-normalize M"'' into B"^' 
by SVD (or QR) to maintain the desired normalization structure. The remaining matrices 
of the SVD are multiplied to the M°''-^ to the left, which will be the starting guess for the 
eigensolver for the next site. Build iteratively the R expression by adding one more site. 
Move on by one site, ^ — > ^ - 1, and repeat. 

• Repeat right and left sweeps, until convergence is achieved. Convergence is achieved if 
energy converges, but the best test is (using MPO) to consider {\p\H^\\p) - {{ip\H\\p)'f- to see 
whether an eigenstate has been reached; this expression should approach as closely as 
possible. 

If we call matrices A, B, M depending on their normalization (M always being the one on the 
site currently attended to), and giving them an subscript index / to label the number of updates 
by the eigensolver they have undergone, the algorithm would formalize as 

MqBqBqBqBqBq 
"^'■^ Ml BoBqBoBoBq A 1 MoBoBqBoBq 
''^ A I M, BoBoBoBo ^-^^ A i A i MqBqBoBo 
''^ A 1 A I Ml BqBqBq ^-^^ a 1 a 1 a 1 MoBoBo 
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AiAiAiAiMiBo AiAiAiAiAiMo 

diag S VD 

-> AiAiAiAiAiMi -> AiAiAiAiMiBi 

diag S VD 

-> A1A1A1A1M2B1 ^ A1A1A1M1B2B1 

t/iflg S VD 

^ A1M2B2B2B2B1 ^ M1B2B2B2B2B1 



and again moving from left to right, starting with a diagonahzation step. 

In this iterative process, the energy can only go down, as we continuously improve by varying 
the parameters. Two problems occur: starting from a random state, the guesses for the M""' in the 
iterative eigensolvers will be very bad in the initial sweeps, leading to large iteration numbers and 
bad performance. Moreover, we cannot guarantee that the global minimum is actually reached 
by this procedure instead of being stuck in a non-global minimum. 

One way of addressing the first issue is to start out with infinite-system DMRG to produce 
an initial guess; an optimal MPS version of infinite-system DMRG is discussed in Section [TOl 
While this initial guess may be far from the true solution, it will usually fare much better than a 
random starting state. Moreover, one can try to balance the number of iterations (high in the first 
sweeps) by starting with small D, converge in that ansatz class, enlarge D and add zeros in the 
new matrix entries, converge again, and so on. When D gets large, the guess states will hopefully 
be so close to the final state that only very few iterations will be needed. It turns out, however, 
that starting with too small D may land us in a non-global minimum that we will not get out of 
upon increasing D. Quite generally, as in DMRG, one should never calculate results for just a 
single D, but increase it in various runs until results converge (they are guaranteed to be exact in 
the D — > 00 limit). 

If we are looking for low-lying excited states instead of a ground state, two typical situations 
occur: (i) The excited state is known to be the ground state of another sector of the Hilbert 
space decomposed according to some good quantum number. Then the calculation is just a 
ground state calculation in that different sector, (ii) The excited state is the first, second, or 
higher excitation in the sector of the ground state. Then we have to calculate these excitations 
iteratively, and orthonormalize the state with respect to the lower-lying states already identified; 
this clearly limits the approach to a few low-lying excitations. The place where the algorithm 
is to be modified is in the iterative eigensolver; e.g. in the Lanczos iterations, the next Lanczos 
state generated is orthonormalized not only with respect to the previous Lanczos states, but also 
already constructed eigenstates of the Hamiltonian. This is a standard extension of the Lanczos 
algorithm. 

The variational MPS algorithm just introduced is quite prone to getting stuck. How this is 
going to happen, actually depends a bit on how initial states are chosen in the procedure: Assume 
that, as is the case for the anisotropic Heisenberg chain, there is a quantum symmetry with some 
commuting operator [H, O] - 0, in this case the total magnetization operator M - Yji^'- , giving 
rise to magnetization M as a good quantum number. Then initial states fall into two categories, 
whether they are eigenstates of M or not. The latter will generally be the case if the state is 
chosen randomly; the former is the case if it is generated by infinite-system DMRG or its MPS 
variant. 

Decomposing the Hilbert space into eigenstates of magnetisation, '7/ = ®m'J^m, we can write 
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initial states as 

|(A> = ^tAM|M) me'HM. (212) 

M 

The ground state we are looking for is \iff()} e 'Hi^; the initial state will then have either arbitrary 
i^M or iJ/M = if M 9^ M (assuming that we don't run into the desaster of offering an initial 
state in the wrong symmetry sector). Let us assume that in the first case, sweeping will eliminate 
contributions from the wrong symmetry sectors; if they don't, the variationally optimal state can 
never be reached anyways because wrong admixtures survive. As an iterative ground state search 
by e.g. Lanczos is an optimized version of the power method lim„^co H"\iJ/) for finding the largest 
eigenvalue and associated eigenstate, one can show that in the full Hilbert space wrong symmetry 
sectors will definitely be projected out. In our algorithm, this iterative projection proceeds in a 
highly constrained state space and might not be as eflicient, as it looks at various wave function 
components sequentially. As random starting states are very inefficient, I cannot report on a lot 
of practical experience here. In any case, once we arrive in a well-defined symmetry sector, we 
will have, for any Schmidt decomposition - Y^ac ^ar\o(}A\ae}B, that each of the states will 
have a good quantum number (superpositions of different quantum numbers lead immediately 
to a contradiction to a global good quantum number), namely and such that + mf - 
M, where I have simplified indices. Taking the mf , for example, they will be distributed over 
some range, say 1 state with magnetization m, 3 states with magnetization m', 5 states with 
magnetization m" and so forth. As I will show next, this distribution stays^xet/in further sweeps. 
This means that if it does not correspond to the distribution that the variationally optimal state 
would yield, it can never reach that state. In the random state approach one may hope that the 
slow elimination of other total magnetizations "eases" us into the right distributions but there is 
no guarantee; in the infinite-system approach one has to hope that this warm-up scheme produces 
the right distribution right away, which is quite unlikely to happen. 

The reason why the distribution stays fixed can be seen from the SVD of M^/ ^ to carry out 
one (for example) left-normalization step: reshaping matrices M°'' into some and applying an 
SVD gives at most D non-vanishing singular values; the right-singular vectors in V ' are nothing 
but the eigenvectors of ^F'^F, which is block-diagonal because the states lai)^ have good quan- 
tum numbers. The right singular vectors (eigenvectors) therefore encode a basis transformation 
within blocks of the same quantum number, hence the number of states with a given quantum 
number remains the same, and so does the number of states with a given quantum number in the 
other part of the system because of the matching of quantum numbers required in the Schmidt 
decomposition. 

Various ways of getting out of this potential trap have been proposed. The first one is to 
modify the algorithm to consider two sites at the same time, just as in conventional (two-site) 
DMRG; we will discuss its MPS implementation in the next section. While this approach is 
slower (roughly by a factor of d), it offers a slightly enlarged ansatz space with a subsequent 
truncation that allows the algorithm to be more robust against the danger of getting stuck in local 
energy minima in ground state searches. In particular, the enlarged ansatz space of the two-site 
algorithm allows a reshuffling of the quantum number distribution due to the truncation. Once 
this is converged, one may switch to the single-site algorithm, as proposed by Takasaki et al. 
ll45h . although it is not at all clear that this leads strictly to the optimal outcome ll65t . 

Much better, there is a procedure by White 1,46] that protects reasonably against trapping and 
ensures reshuffiing. It is crucial for a reliable single-site DMRG (or variational MPS, which we 
will show to be identical) algorithm and turns it into the state of the art form of the method. It 
starts from the observation that quantum numbers of a subsystem A are changed by quantum 

69 



? T T T 




Figure 43: (H^'li//))"'-' represented graphically: with the exception of one W-tensor and one 'F-tensor, all the 



contractions have already been computed to obtain L. 



fluctuations due to those parts of the Hamihonian that connect A to the rest of the system. We 
therefore have to consider the structure of H\i//} in more detail. 

Consider i//aii.ar oft^f energy optimization in the single-site algorithm. We can also write 



br 



B 



(213) 



where 



(214) 



cr.cr'eA. 



cr,cr'eB 



such that there are Dw terms in this sum. If we think in terms of block states, we would like to 
know which new states can be reached on A» by the action of i/^V Projecting the result of this 
action onto the A»B basis it will read 



b\ I, a. 



/ i «f-2,af-l bc-2,he-l flf_2.«f-i 



bt-i,be a'i_i,ar 



(215) 



which is just 



(216) 



usmg lI"'"'-' fromEq. ( fT92] i. as can be seen graphically in Fig.|43] This indicates that the actual 
cost of computation is very low, because we have already done the most complicated part. 

Now we would like to include the states generated by i/^* into the search for a good basis 
for A». Here, DMRG offers the possibility of multiple-state targetting. The conventional algo- 
rithm would now proceed by calculating p^' = TiB\i//){4r\ or P(af_,o-e),(ai,c-'^) = "^aLuat^'il'^ , ' 
finding the eigenvalues (squares of the singular values), the eigenvectors (left singular vectors), 
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truncation and so on. But we can look at a modified density matrix, which takes also into account 
the new terms as 

p^- = TrBm4f\ +aj] TrgHt- mMH^', (217) 

bi 

where a is a small number, say lO"'*, giving a little weight to contributions that the conventional 
algorithm may miss. The price paid is that at the end of the spectrum a few very-small weight 
states from will drop out. Upon multiple sweeping, a will be taken slowly to zero. 

The new density matrix is diagonaUzed, and truncation according to the D largest eigenvalues 
is carried out, yielding a {dD X D) matrix of orthonormal columns, U{a,_fa,),a, ^ ^fl/_,,af . and we 
continue on the next site; as we are not using the eigenvalues of the modified density matrix 
beyond their relative ordering, it does not matter that they do not sum up to 1. For predicting the 
next for the large sparse matrix solver, we use the DMRG prediction formula derived in the 
next section, 

vi/o-f+i _ A"''^ xu^'e nO'e+i 

aeaui a,_ia^ a^,af+i' K^'^°J 

ae-ia-ea'f 

Otherwise, everything remains the same. The additional numerical cost and programming effort 
is minimal for an algorithm that often converges much faster and is much less prone to getting 
stuck at a non-optimal result. 

6.4. Conventional DMRG in MPS language: the subtle differences 

How can the previous approach be related to conventional DMRG? The essential answer is 
that the MPS approach is identical to finite-size DMRG for OBC, albeit only if we shift one 
site instead of two, i.e. consider "single-site" instead of "two-site" DMRG, where we consider a 
block-site-block configuration A»B instead of a block-site-site-block configuration A»»B. 

Let us first remind ourselves of the key steps of the algorithms, assuming that we are sweep- 
ing to the right: (i) given some configuration A»B (or corresponding configuration AAAMBB) 
and a Hamiltonian H, the ground state is found by a large sparse matrix eigensolver looking for 
the optimal i//ar_itT[a[ (in DMRG) or M'^'^_^ at (in MPS) respectively; analogously for A»»B. (ii) 
Given the ground state, MPS derives a set of left-normalized A-matrices, whereas DMRG finds 
new block states whose structure can be encoded by left-normalized A-matrices. (iii) All algo- 
rithms switch to a new A»B, A»»B ox AAAAMB configuration, where the active center is shifted 
by one site to the right and provide an initial guess for the calculation of the next ground state, 
taking us back to step (i). 

Step (i): Results must be identical if we use the same state configuration and the same Hamil- 
tonian. As DMRG grows the blocks A and B from left and right, and as each block growth step 
A» ^A can be encoded by A-matrices and similarly "B — > B, we conclude that all matrices on 
A are left-normalized and those on B right-normalized, hence the two-site DMRG state takes the 
form 

\^) = Yj \ae-i)A\'re)\(Te^{)\ae^i)B = 2 ' ' . A'^'-^P^B'^'- . . . B^'\cr), 

of-iffff+iflf+i cr 

(219) 

with the obvious change for a single-site DMRG state, *P°"<'^<+i — > *P°"f . This is in perfect agree- 
ment with the mixed-canonical states of the variational MPS approach and we are looking at the 
same state structure. 
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It remains to show that the Hamiltonians are identical, too. Strictly speaking, this is not the 
case: The MPO representation of H we just used is clearly exact. On the other hand, the repre- 
sentation of H in DMRG contains a series of reduced basis transformations, hence is inherently 
inexact. So, the two representations seem unrelated, with an advantage on the MPO side because 
it is exact. But a more careful analysis reveals that on the level of calculating expectation val- 
ues {tfr\H\il/} as they appear in MPS and DMRG ground state searches both representations give 
identical results (they are not identical for higher moments, such as {t//\H^\i//}, where the MPO 
representation is demonstrably more accurate at a numerical cost, see below). 

Both the DMRG and the MPO Hamiltonian contain all terms of the exact Hamiltonian. As we 
have already seen in the application of a Hamiltonian MPO to a mixed canonical state (Sec. l6.2l l. 
the evaluation of the L and /^-objects appearing in the large sparse eigenproblem Eq. (1209b is 
nothing but the sequence of reduced basis transformations occuring in DMRG up to the cur- 
rent A»B configuration. Hence, for {iff\H\4r) (but in general only for this!), both approaches are 
identical. 

Moreover, the calculation H\iff} appearing in the eigenproblem does not have a worse opera- 
tional count than in the corresponding DMRG procedure. To see this, let us focus on our example 
MPO for an anisotropic nearest-neighbour Heisenberg chain. There seems to be a difference in 
efficiency when we consider the double sum over b(-\,b{. From the structure of the W-matrix 
it is clear that for most of the (in the example 25) entries we find zeros, such that we can 
strongly restrict the sum. But this would still give the following count: setting the field 0, for the 
Heisenberg Hamiltonian there are 8 contributions in the DMRG setup: one each for Ha and Hb, 
the parts of the Hamiltonian that act strictly on A and B, and three per block for the three operator 
combinations linking a block and the site. All of them are diagonal in the other block, so there 
are altogether 8 operations of cost O(D^). In the MPO calculation, the double matrix-matrix mul- 
tiplication would naively suggest 16 operations of cost O(D^), for the 8 non-vanishing entries of 
W°''''^r. But then we can exploit the following rules: if bf^i = dw, then there are no operations 
to the left, L^J^^''"" - 6a[_^^a'^ ^, and one operation drops out. Similarly, \fb[ - 1, then there are 

no operations to the right and R'^^'"' - Sa[,a[, and again one operation drops out. Looking at the 
structure of W, all non-vanishing entries meet one or the other condition, and the count halves 
down to 8 operations. Only longer-ranged interactions do not fit this picture, but they would be 
of cost 20{D^) in DMRG as well. 

Step (ii): After energy minimization, variational MPS and DMRG produce (identical) M'^' 
and , or >I'°"f°"f+i. Both methods now seem to proceed differently with the result, but in fact 
do the same: in variational MPS one just shifts one site to the right after an SVD to ensure 
left-normalization, to continue minimizing on the next site. In DMRG one previously carries 
out a density matrix analysis to determine a new (truncated) block basis. But if one carries out 
the corresponding SVD, the number of non-zero singular values (hence non-zero density matrix 
eigenvalues) is limited by D, the matrix dimension: 

mm(dD,d)=D 

ni„«. ^ ^(«.-,<..).«. = Yu KLk^kk{y%ar (220) 

k=\ 

Hence, no truncation happens, and we are just doing a unitary transformation to obtain orthonor- 
mal states for the new larger block A (which is just the left-normalization in the MPS because of 
the link between SVD and density matrix diagonalization). Both formalisms act identically; as 
no truncation occurs, thin QR would do, too. 
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On the other hand, in two-site DMRG the same step reads 



mm(dD,dD)=dD 



(a(-io-f),(o-f+ia(+i) — 



2 



flf+l)- 



(221) 



But we can only keep D states in the new block, hence truncation has to occur! Here is the only 
difference between variational MPS and single-site DMRG on the one and two-site DMRG on 
the other hand. 

Step (iii): In DMRG, after completion of one iteration, the free site(s) are shifted by one, 
leading to block growth of A and shrinkage of B. Here, all methods agree again: in variational 
MPS, the shrinkage of B is simply reflected in the states being formed from a string of B-matrices 
where the leftmost one has dropped off. The growth of A is given by a similar string, where one 
A-matrix has been added. The matrix on the free sites is to be determined in all approaches, so 
nothing is to be said about its normalization. 

Minimization of ground state energy is, as we have seen, a costly large sparse matrix problem. 
As the methods are iterative, a good initial guess is desirable. DMRG has provided some "state 
prediction" for that [lY]. In fact, it turns out that the result of the prediction is just what one 
gets naturally in variational MPS language without the intellectual effort involved to find state 
prediction. 

Let us assume that for single-site DMRG we just optimized deriving a new A"''. Then 
in MPS language the next ^I"^-^*' = SV^B°''*\ where S and V' are from the SVD. In DMRG 
language, we take 

1-^)= X "^aLJ^^-^yAlcriM)^ (222) 

and insert twice approximate identities / = 2 |flf)^ ^(flfl and / = 2 \o'e+\}\ae+i}B B{ae+\\{crc+i\. 
Expressing the matrix elements by A and B matrices, the state now reads 



Z Z 



(223) 



So the prediction reads 



vr/o" (+1 



"f )"f-i ac-\ a^ae-\-\ 



(224) 



But fliis is exactly the MPS ansatz for the next eigenproblem, as A'^'^'^°''B°''*' = A°'''^A°''S V'^B"''*' . 
But this is just S y ' because in the ansatz, cr( is summed over and left-normalization holds. 
Two-site DMRG proceeds by analogy and is left as an exercise for the reader 

While this clarifies the relationship between variational MPS, single-site DMRG (the same) 
and two-site DMRG (different), it is important to note that the different ways of storing infor- 
mation more implicitly or more explicitly implies differences even if the algorithms are strictly 
speaking identical - the fact that in one formulation prediction is trivial and in the other is not 
already gave us an example. But there is more. 

(i) In DMRG, the effective bases for representing the states and the Hamiltonian or other 
operators are tied up. This is why concepts such as targetting multiple states arise, if we consider 
several different states like the ground state and the first excited state at the same time. One then 
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Figure 44: "Exact" calculation of the expectation value of H^ : the Hamiltonian MPO is repeated twice and sandwiched 
between \\p) at the bottom and at the top. 

considers mixed reduced density operators 

Pa = Yj ai^^BWi){n (225) 

i 

with li/r,) the target states and < a,- < 1, 2, ff,- = 1, to give a joint set of bases for all states 
of interest. This can of course only be done at a certain loss of accuracy for given numerical 
resources and for a few states only. At the price of calculating the contractions anew for each 
state, in the MPO/MPS formulation, the state bases are only tied up at the level of the exact full 
basis. MPO/MPS formulations therefore acquire their full potential versus conventional DMRG 
language once multiple states get involved. 

(ii) Another instance where the MPO/MPS formulation is superior, albeit at elevated numer- 
ical cost, is the calculation of the expression {i//\H^\i//}, which is interesting e.g. in the context 
of estimating how accurately a ground state has been obtained. In the MPO formalism, it can 
be done exactly up to the inherent approximations to Itfr) by contracting the network shown in 
Fig. |44] It would of course be most economical for the programmer to calculate H\tfr} and take 
the norm, two operations which at this stage he has at hand. The operational cost of this would be 
O(LD^D^d^) for the action of the MPO and 0{LD^D\fd) for the norm calculation. The latter is 
very costly, hence it is more efficient to do an iterative construction as done for {il/\H\4r}. Let me 
make the important remark that dimension is only the worst case for H^|74]: Writing out the 
square and introducing rules for the expression leads to more efficient MPOs, whose optimality 
can be checked numerically by doing an SVD compression and looking for singular values that 
are zero. Our anisotropic Heisenberg Hamiltonian takes Dw - 9 instead of 25 for H^. For higher 
powers, the gains are even more impressive, and can be obtained numerically by compressing an 
explicit MPO for H" with discarding only zeros among the singular values. 

In a DMRG calculation, there would be a sequence H{H\i(f}), in the DMRG block-site basis 
as shown in Fig.|45] The point is that before the second application of H, a projection onto the 
reduced block bases happens, which is not the identity and loses information. 

What does the comparison MPS and DMRG imply algorithmically? First of aU, the trunca- 
tion error of conventional DMRG, which has emerged as a highly reliable tool for gauging the 
quality of results, is nothing but an artefact of the somewhat anomalous two-site setup. In varia- 
tional MPS or single-site DMRG it has to be replaced by some other criterion, like the variance 
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Figure 45: DMRG calculation of the expectation value of H-: the Hamiltonian MPO is applied once to \if/) in step (1) in 
the DMRG basis, i.e. the result is projected onto the reduced bases, yielding some <I)°"f . This in turn replaces 'F""' in the 
second application of i? in step (2). Ultimately, the result is projected on the block bases. 

of the energy. Second, while all the approaches are variational in the sense that they are looking 
for the lowest energy that can be achieved in a given type of ansatz, it varies from site to site 
in two-site DMRG (because of the 'P°'°' anomaly in the ansatz), the ansatz stays the same all 
the time in single-site DMRG, which is conceptually nicer That this comes at the expense of 
potential trapping serves as a reminder that the mathematically most beautiful does not have to 
be the most practical. 

7. Time evolutions (real and imaginary) with MPS 

The calculation of the action of operators like e '''' or e '^'^ on quantum states is of central 
interest in quantum mechanics, for real-time evolutions of quantum states and for quantum sta- 
tistical mechanics; /3 can be interpreted as an imaginary time. It is one of the most attractive 
features of MPS that such real or imaginary time evolutions can be encoded very neatly and 
efficiently. This holds both for pure and mixed states, important at finite temperature. In the 
following, I will focus on time evolution based on a Trotter decomposition of the evolution op- 
erators 1 47, 48, 4^ S^, 51], explaining first the Trotter decomposition and the structure of the 



algorithm for pure states, then the representation of the Trotter decomposition by MPOs. After 
this, I will discuss the changes necessary for the simulation of the dynamics of mixed states. 

7.1. Conventional time evolution: pure states 
7.1.1. Trotter decompositions of time evolution 

Let us assume that H consists of nearest-neighbour interactions only, i.e. H - 2, /z,, where 
hi contains the interaction between sites / and / -i- 1 . We can then discretize time as f = Nt with 
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T — > 0, A'^ — > 00 and (in the most naive approach) do a first-order Trotter decomposition as 



er^-" = e-iAirg-iA2Tg-iA3T _ _ _ g-iAt_3Tg-ife_2Tg-iftz._,r ^ (226) 

which contains an error due to the noncommutativity of bond Hamiltonians, [hi,hi+i] in 
general; higher order decompositions will be discussed later in this section. All time evolutions 
on odd (e"''^"""^) and even (e"'^"™-'') bonds respectively coimnute among each other, and can be 
carried out at the same time. So we are looking for an MPO doing an infinitesimal time step on 
the odd bonds and for another MPO doing the same on the even bonds. 

As any operator is guaranteed to be MPO-representable, let us assume for a moment that in- 
deed we can construct these representation of infinitesimal time steps efficiently (see next section 
for the explicit construction). As we will see, the maximum bond dimension of the infinitesimal 
time step MPOs is cfi because the dimension of e'^'^ is {d^ x d^). The apphcation of the in- 
finitesimal time step MPOs thus increases the bond dimensions from D up to d^D. Repeated 
applications of the infinitesimal time evolution MPOs leads to an exponential growth of the ma- 
trix dimensions, which therefore have to be truncated after time steps. 

The resulting time evolution algorithm takes a very simple form: starting from \tj/(t = 0)), 
repeat the following steps: 

• Apply the MPO of the odd bonds to \tfr(t)). 

• Apply the MPO of the even bonds to e~'"°^''\ifr{t)). 

• Compress the MPS \t//(t + r)) = e"'^"="^e~'^°'"^|^(?)) from dimensions d^D to D, monitor- 
ing the error Obviously, one may also allow for some compression error (state distance) 
e and choose a time-dependent D: it will typically grow strongly with time, hmiting the 
reachable timescale. By analogy to the ground state calculations, all results should be 
extrapolated inD— >ooore— >0. 

After each time step, we may evaluate observables in the standard way, (0(?)) = {tl/it)\0\t//{t)). 
But we can do more: we can calculate time-dependent correlators as 

{dit)P) = {,fr\c^'"'dc-^'P\0 = (miom)) (227) 

where Mt)) = e-^'\0 and Mt)) = e-''"'P\il/). If we take e.g. O = and P = S% we can 
calculate {S^{t)S^^ and by a double Fourier transformation the structure function 

S''(k,oj)^ fdtY,(S'i(t)Sl„W'"e-"'', (228) 

where I have assumed translational invariance and infinite extent of the lattice for simplicity of 
the formula. 

A simple improvement on the algorithm given above is to do a second-order Trotter decom- 
position 

g-iHr ^ g-iH<Kjdr/2g-iHevenTg-iHoddr/2 _^_ 0(t^), (229) 

where the error per timestep is reduced by another order of t. If we do not do evaluations after 
each time step, we can group half steps, and work at no additional expense compared to a first- 
order Trotter decomposition. 
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Figure 46: A Trotter step: On all odd bonds, an (infinitesimal) bond time evolution is carried out. This merges two sites, 
such that the simple product form of MPS is lost at first sight, but the time evolution can be translated into MPOs. As the 
time evolution factorizes, the MPOs have dimension 1 on all even bonds (thin lines). 



A very popular implementation of a fourth order- Trotter decomposition that origi nates in 
quantum Monte Carlo would be given by the following formula due to Suzuki lll02[ll03n : 



,-iHT 



t/(Ti)t/(T2)t/(T3)f/(T2)t/(Ti), 



where 



and 



U(t) — e"'^"'''''^'^^e '^""'"'e^'^"'"'''"'^^ 
1 



4-41/3 



2ti - 2t2 



(230) 
(231) 
(232) 



Even smaller errors can be achieved at similar cost using less symmetric formulae 111 0411 . This 
completes the exposition of the algorithm (an error analysis will be given after the other methods 
for time evolution have been explained), and we now have to construct the MPOs. 

7.1.2. MPO for pure state evolution 

Let us consider the Trotter step for all odd bonds of a chain: 



-ih\T , 



)e"''''-'^|(A); 



each bond-evolution operator like e ''^ takes the form 2 



t' ^ 



(233) 

|o-iCr2)<cr'jCr^|. Both 



in the pictorial and the explicit mathematical representation it is obvious that this operator de- 
stroys the MPS form (Fig.l46]l. 

It would therefore be desirable to have (9°"i°"2 °"'i'^2 in some form containing tensor products 
(9°"''°"'i (gi (9°"^'°"2, to maintain the MPS form. To this purpose, we carry out the procedure for 
decomposing an arbitrary state into an MPS, adapted to an operator (two indices per site). It 
works because there are so few indices. One reorders O to group local indices and carries out a 
singular value decomposition: 



O' 



(iTio-'),(a-20-;) 



k,{o-20-') 
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Figure 47: A complete first-order Trotter time step (odd and even bonds). Fat and tliin lines correspond to dimensionl 
and > 1 on MPO bonds. The MPOs in the top line on the first and last site ai'e trivial scalar identities 1. 




where f/^ ' ' = U(a-xcr\ ),k sjS k.k and f/i' ' = ylS^k(V^)kXcr2a-',)- In the very last step of the deriva- 
tion, we have introduced a dummy index taking value 1 to arrive at the form of an MPO matrix. 
The index k may run up to d^, giving the bond dimension Dw of the MPO. 

The MPO representing the operator in Eq. ( |233] |. U"''"'''ir'°'^ V^-'^'^lT'''^' . . ., factorizes on 
every second bond, as do the original unitaries. If one site does not participate in any bond 
evolution, we simply assign it the identity unitary as a (1 x l)-matrix: f^f = Sa-,cr'- Then 
the global MPO can be formed trivially from local MPOs. The MPO for time evolution on all 
odd bonds would read UUUUUU . . ., whereas the even-bond time step reads lUUUUUU . . . I 
(Fig. 113. 

7.2. Conventional time evolution: mixed states 
7.2.1. Purification of mixed states 

Finite temperature calculations can be carried out based on the purification of an arbitrary 



mixed quantum state 115 lH : if we consider a mixed state in physical space P formed from or- 
thonormal states, we can interpret it as the result of a partial trace over a Schmidt decomposition 
of a pure state on PQ, where Q is an auxiliary space: 

r r 

= ^ sl\a)p{a\p ^ = ^ Sa\a)p\a)Q pp = TrQ|i/r)<(/r|. (234) 

n=l fl=l 

The auxiliary state space can simply be taken as a copy of the original one, so finite-temperature 
density operators on a chain can be expressed as pure states on a ladder (see Fig.l48ll. 
To calculate a thermal density operator /jg - Z(J3y^e^^", Z(J3) - Tr^e"^^, we write 

p/s = ZiPT'c-l"" = Z(/?)-'e-^«/2 . ; . g-M/2_ (235) 

The identity / is nothing but Z(0)po, the infinite temperature density operator times the infinite 
temperature partition function. Assume we know the purification of po as an MPS, |iA/3=o)- Then 



(Z(Q)IZ{l3))t-P"l^ ■ Trel^oX^ol ■ ^'^"^^ = {Z{Q)IZ{J3))i:rQt-P"^^mmt-P"l\ (236) 
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Figure 48: Schematic representation of finite-temperature simulations: instead of a chain, one sets up a ladder with an 
identical copy of the chain. Physical sites have odd. auxiliary sites have even labels. Equivalent sites on the physical and 
auxiliary leg are linked by maximally entangled states. To reach inverse temperature j8, an imaginary time evolution is 
carried out up to "time" -i/?/2. 



The trace over Q can be pulled out as the Hamiltonian does not act on Q. But the result means 
that we have to do an imaginary time evolution 

l<A/3> = e-^'^/Vo). (237) 

Expectation values are given by 

(O)^ = TxpOpp = {Z{0)IZ{pmpdTxQ\ilfp){il,p\ = {Z{0)IZmWO\ilfp). (238) 

(Z{0)/Z{fi)) = {d^/Z{/3)) may seem difficult to obtain, but follows trivially from the expectation 
value of the identity, 

1 = </)^ = Tvppf, = (Z(0)/Z(/?))TrpTre|^;,)<^^| = (Z(0)/Z(;8))<^^|^^), (239) 
hence Z(J3)/Z(0) - {ij/^lilf^}, or, in complete agreement with standard quantum mechanics. 

But this takes us right back to expressions we know how to calculate. All we have to do is to find 
\i(f[)}, carry out imaginary time evolution up to -i/?/2, and calculate expectation values as for a 
pure state. We can even subject the purified state to subsequent real time evolutions, to treat 
time dependence at finite T. 

We can also do thermodynamics quite simply, as Z(J3)/Z{Q) is given by the square of the 
norm of and Z(0) - d^. This means that we can obtain Z(J3) by keeping the purified state 
normalized at all temperatures and by accumulating normalization factors as temperature goes 
down and /3 increases. From Z(J3), we have F(J3) - -y6"' lnZ(J3). At the same time, U(j3) - 
{H)l3 - {t///j\H\iff/s}. But this in turn gives ns S(Ji) - fHUifi) - F(J3)). Further thermodynamic 
quantities follow similarly. 

The purification of the infinite temperature mixed state is a simple MPS of dimension 1, 
because it factorizes (if we take one ladder rung as a big site): 

1 . /I .f^ 

P0 = ^/ = (^/) . (241) 
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Figure 49: The time-evolution of mixed states on a chain can be seen as that of a pure state on a ladder, where the physical 
sites sit on the first leg and additional auxiliary sites on the second leg. This ladder is mapped to a chain. As the time 
evolution acts only on the physical states, next-nearest neighbour interactions arise. 



As po factorizes, we can now purify the local mixed state on each physical site as a pure state on 
rung / , to get |i/',o), then |i^o) - |iAi,o)|i/'2,o)|iA3,o) • ■ •, a product state or an MPS of dimension 1. If 
we consider some rung / of the ladder, with states \o-}p and |cr)Q on the physical site 2/ - 1 and 
the auxiliary site 2i, we can purify as follows: 



= ^ i p\cr){o-\p = Ttq 



1 V 1 



(242) 



Hence the purification is given by a maximally entangled state (entanglement entropy is log2 d), 

l<A,o> = J] ^k>pk>e. (243) 

It is easy to see that one can carry out local unitary transformations on both P and Q separately 
that leave that structure invariant. For example, for the purification of a spin- 1/2 chain it is 
advantageous to use the singlet state as local purification. 



l-A/.o) = ^[| IpIq) - 1 ipTe)], 

V2 



(244) 



in case the program knows how to exploit good quantum numbers: this state would allow to 
conserve total 5=0 and 5^ = at the same time. In this case, the four A-matrices would read 



(245) 



and the purified starting state |i/'o) for /3 - is now given by a product of singlet bonds on a 
ladder. In fact, a S VD of reshaped matrix Aa-^,j> allows us to introduce truly site-local A-matrices, 
which have dimension (1 x 2) on odd and (2 x 1) on even sites: 



A^-'-' = [1 0] A^='-' = [0 - 1] A''-' = [0 1/V2]^ A^-< = [1/ V2 0]^ . 
In order to apply the pure state time evolution algorithm, it remains to find the MPO. 



(246) 
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7.2.2. MPO for mixed state evolution 

The ladder appearing in mixed state simulations can be mapped to a chain (Fig.l49b. where 
the physical Hamiltonian acts only on the odd sites, 1, 3, 5, . . ., and the auxiliary sites are even, 
2,4, 6, Then the only non-trivial time-evolution connects (1,3), (3, 5), (5, 7). There are sev- 
eral ways of dealing with such longer-ranged interactions, one explicitly constructing the longer- 
ranged interaction, the other using so-called swap gates, reducing it to a nearest-neighbour inter- 
action. 

The direct MPO description of the "longer-ranged" interaction (1,3) involves necessarily a 
non-trivial tensor on site 2, whereas the site 4 is inert. Similarly for (3, 5), there is a non-trivial 
tensor on 4, but site 6 is inert. This suggests a Trotter decomposition (1,2, 3, 4), (5, 6, 7, 8), . . . in 
the "odd" and (3, 4, 5, 6), (7, 8, 9, 10), ... in the "even" steps. 

The four-site evolution operator on sites 1 through 4 then reads 

^(iTio-20-3iT4),(o-',a-;a-',o-p _ Q(a-iU2Cr3),(cr\a'-^a-'^) ^ (247) 

and we can build the three-site unitary with only a slight modification of the two-site unitary 
which contains the actual physical time evolution: 

^(o-io-2iT3),(o-;o-;o-',) _ Q(a-iai),(a-[cr',) ^ ^ (248) 

This three-site unitary is now subjected to two SVDs. For the notation, we first shift down the 
indices and reorder sitewise. Reshaping with subsequent SVDs then iteratively isolates 0-3,0-^, 
0-2,0-'^, and <j\,a-\: 

— -P(o-iit',0-20-2),(o-30-3) 

~ ^ ^(o"lO"',a"2a"2).*2'^L'i:,(^23-'*2.(o'30'3) 
k2 

~ XI ^('^i'^l)''^i'^lI!*:/^12^*i'(°"2<*2)'^ij^A;2^^23^*2,(tr30-;,) 
k, ,k2 

Zj ^^l,ki ^^hM ^^k2,l 



where, with the introduction of dummy indices and the inert tensor on site 4: 



Vc^(v;2).„(.2<*2)V^ (250) 

" V^^(^23)*2,(a-3tr;,) (251) 



kiM 



= ^-4.< (252) 

From this, MPOs for the entire chain can be formed as for the pure state time evolution. We 
have done nothing but the iterative decomposition of an MPO on 4 sites. Again, this is still 
manageable, as only 4 sites are involved. 
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Obviously, it is a straightforward step to write an evolution operator acting on all four bonds 
(1, 2), (1, 3), (2,4) and (3,4) and subject it to a similar sequence of SVDs, which would allow to 
consider pure state time-evolution on a real ladder. 

An alternative approach to carry out interactions beyond immediate neighbours is provided 
by the use of swap gg fe^ lilOOi] . Let us take the example of a real ladder with interactions on 
four bonds, two of which [(1, 3) and (2,4)] are next-nearest-neighbour interactions. But if we 
swapped states on sites 2 <-> 3, they would be nearest-neighbour interaction. Time-evolution on 
the ladder would then be done as follows: (i) evolve bonds (1,2) and (3,4); (ii) swap states on 
sites 2 and 3, (iii) evolve "new" bonds (1,2) and (3,4) with the evolution operators for "old" 
bonds (1,3) and (2,4) and (iv) swap states on sites 2 and 3 once again. In other situations, more 
astute schemes need to be found, preferrably generating a sequence of swaps between nearest 
neighbours. The swap operator for sites i and j is simply given by 

Sij^ Z ^"'"^"■"V;cr^)<(r;cr;| S'^^''^'^'^' ^ 6^,^^6^^,^^, (253) 

o'iCr'.a'jO''. 

is unitary and its own inverse. It swaps the physical indices of two sites in an MPS; for swaps 
between nearest neighbours it is easy to restore the original form of the MPS: assume that the 
MPS is left-normalized; a unitary applied to sites / and /+ 1 affects this only on these two sites. In 
particular, the orthonormality of block states |flj:<,)A and |a,+i)A is not affected. If we introduce a 
matrix M(a,_,^,^,),(^,,a,^,) = Za, '4|,'1'[;+;,A[,':|„'J^^', we can form M(a,_ta-:Ua-MaM) = ^,,_^,),(<r,,a,^,) and 
carry out an SVD, where U(ai_i(ri),ai yields a new left-normalized Agi_^ g. and S ai,ai{V^)ai,(a-Mai+i) a 
new left-normalized Ag'*^.^^ . That the latter is left-normalized follows from the left-normalization 
of A""' and the maintained orthonormality of the |fl,+i)yi. 

Let me conclude this outlook on beyond-nearest-neighbour interactions with the remark that 
using MPO allows also other Trotter decompositions, e.g. decomposing the Heisenberg Hamil- 
tonian in its x, y and z-dependent parts, useful for long-range interactions i72ll . 

7.3. tDMRG and TEBD compared to MPS time evolution: the little differences 

A bit before time evolution with MPS (tMPS) was developed, two other algorithms were 
introduced to simulate the real-time dynamics of one-dimensional quantum chains, time-evolving 
block decimation (TEBD) Hill and real-time or time-dependent DMRG (tDMRG) HIs^]. 
Both algorithms are also based on MPS, but are different from tMPS, when one looks more 
closely. Before I get into that, let me stress however that all of them are based on the idea of 
time-evolving an MPS which was first put forward in i47l l48ll and therefore are minor variations 
on a theme. tDMRG and TEBD are mathematically equivalent, i.e. should for exact arithmetic 
give the same results, whereas numerically they are clearly distinct algorithms, both carrying out 
operations that have no counterpart in the other method, with their respective advantages and 
disadvantages. Let us discuss first tDMRG, because its language is closer to that of tMPS, and 
then TEBD, to see how important (or unimportant) the little differences are. 

7.3.1. Time-dependent DMRG (tDMRG) 

The decomposition of a global time-evolution on an entire lattice into a Trotter sequence of 
infinitesimal time-evolutions on bonds is the same for all three algorithms discussed here. Let us 
therefore focus on one infinitesimal time-evolution e"''''*''" on sites £ +1 and ( + 2. The evolution 
operator expressed in the local basis is given by 
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The current state if/ is given in the two- site DMRG notation with left- and right normalized 

matrices as 

\if/) = Y^A°'' ... A'^'^V'^-^'^'^^B^'*^ . . . B°''-\cr). (255) 
cr 

The time-evolution turns >i"^'+i'^f+2 into 

<xr= 2 f/<"'^'"^*^^'("''-"'-^<ri!i". (256) 

This, together with the A and B-matrices defines a valid DMRG state we call In order to 
make progress, namely to apply eT^'*^'^ on the next pair of sites, we have to bring the state into 
the form 

10) = ^ A°"> . . . A'"*^^'^'*^'^'*'B^'*' . . . B^'-\(r). (257) 

(T 

The changes can only concern sites £+1 through ^-1-4: on the first two sites because of the action 
of the evolution operator, on the last two sites because they are brought into DMRG form. Let us 
first generate the new A-matrices on sites (+ \ and 1 + 2: We reshape ^Z'^aZi' ~ ^(nfff+i) (o-f+2af+2) 
and subject it to an SVD (DMRG traditionally does this by a density matrix analysis and the 
DMRG prediction when shifting sites, leading to the same result): 

(258) 

oe+i 

U can immediately be reshaped into a valid A-matrix, but has column dimension up to dD, 
which has to be truncated down to D while maintaining the best approximation to the MPS of 
dimension dD. The answer is provided as always by keeping just the D largest singular values 
and shrinking the matrices U, S , V ' accordingly. Here lies the approximation of the method 
(beyond the obvious Trotter error). This done, we reshape as 



2' 



^a^Me^,^ aui,a(+i (^^)a(+1',a<+2 (259) 



and form 3) shifted by one site as 

fl,a-f+20-,+3 _ V c fT/t^o-m DO-f+s f2601 

ar+2 

But we have to shift by another site, which we achieve by reshaping Oa/*, ^^^l^ as ^(a(+i(Tr+2),{o'r+^a(+i)^ 
carry out an SVD as done before, keep the states corresponding to the D largest out of dD singular 
values, reshape, note down A'^'*^ and form 3) shifted by two sites as 

^af+2,ac+4 / J " af+2 . Of +2 V*^ /flM-Of+s af+3.af+4' K^^'^) 

The second SVD and the associated truncation down to D singular values does not lose further 
information, because there are at most D non-zero singular values, although formally there could 
be dD of them. The reason is that before the time evolution on sites i + I and { + 2, the Schmidt 
rank across the bond { + 2 was at most D (due to the MPS construction). The Schmidt rank of 
two states is however identical if they are related by a unitary transformation that acts on either 
part A or part B. But the infinitesimal time-evolution was a unitary on part A. 

We can now continue with the next infinitesimal local time-evolution step, in the spirit of 
tMPS. 
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73.2. Time-evolving block decimation (TEBD) 

Here, we assume that we have in the FA-notation, 

(T 

This state can be immediately connected to the two-site DMRG notation. In particular, 

\j/o"f+io"f+2 — f^i]Y°''*^ K^^^^'^T'^'*-AS''^^^ (263) 
This is identical to the DMRG ^, so is the evolution operator {/(°"f+i°"'+2)'(°"f+i°"f+2), hence also 

Z f/^"^^'"'^^^'<"->"'-^^<iCr- (264) 

In order to proceed, the FA-notation has to be restored on the two active sites. In perfect analogy 
to tDMRG, one obtains by SVD 

*^(flfO-f+l)>(o'f+2af+2) ~ Z ^(afD'f+l)>af+l^af+i,af+i(^'^)af+l>(o'f+2af+2)' (265) 

What is missing are A^'^' and A'^^^^ . We therefore write (reshaping U and V ' and omitting the 
a-indices) 

Now, as in tDMRG, there are up to dD singular values in At'^^'', which we truncate down to the 
D largest ones, just as in tDMRG, also truncating the neighbouring matrices accordingly. We 
now introduce 

F°"';' = (At^l):' m'*' V:'*\ = V°:'*^l (At^+21):' „ (267) 

and obtain 

(j)0-f+ia-f+2 _ AMr'^''+' a'^"'"'-'F°"'*- A'-''^^' (268) 

back to the canonical form. In order to consider the time-evolution on the next bond, we have to 
carry out no SVDs, but just group 

vpa"f+30-(+4 _ a[''+^1f""^+-' A'^"'""'-'F°''+''A'-''^^' (269) 

and continue. As in tDMRG, no loss of information is associated with this step, but this is more 
explicit here. 

When D becomes very large in high-precision calculations, singular values will tend to be 
very small, and dividing by them is, as mentioned previously, a source of numerical instability. 
In the context of the thermodynamic limit iTEBD method, which we will dis cuss later, Hastings 
has proposed an elegant workaround that comes at very low numerical costflOS'], but it can be 
easily adapted to finite-system TEBD. Let us assume that we start with a state in representation 
(12621 ). We then group all pairs into right-normahzed Z?-matrices, 

B'^' = F'^'AI'1, (270) 

but remember the A^'^ for later use. We then form 

^f+io-ui _ jjCTi^i — /^[t+UY^l+lJi^V+^l (271) 
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hence 'P'^'+i^'+2 ^ j^ie]^!' 



|7r70"f+io-(+2 



. We carry out the time-evolution on *P 



to obtain 




;o"f+lO"(+2),(cr'f+]0"f+2)>J? 




(272) 



Then 0^^+'^^+^ = A^^l^' 



. As before, we carry out an SVD on (l)°'f+i°"f+2^ to obtain 



Truncating down to the D largest singular values, we have found the new A^^^'', to be retained 
for further usage, and the new B"''+-. The new is given by 



hence costs a simple matrix multiplication; divisions have been avoided. For the last equa- 
tion, we use right-normahzation of hence 2^,^^ <I)'^'*"^'+^B°"'+="'' = A'^^*' A'^+'J. At the 
same time, B°''*' = y'^m^V+U ^ (A[''J)"'A'^f^' A^^+'l. Combining these two identities with 
(j)0-f+ifr,+2 - /^inQ) - gives the result. 

7.3.3. Comparing the algorithms 

Comparing TEBD and tDMRG step by step, one sees immediately the complete mathemat- 
ical equivalence of the methods. The second SVD in tDMRG does nothing but shifting the 
boundary between left- and right-normalized matrices, which in TEBD is simply achieved by 
a rebracketing of F and A. Nevertheless, there are differences: tDMRG carries out two costly 
SVD decompositions (or density matrix analyses, which is equivalent) per bond evolution, where 
TEBD does only one. On the other hand, TEBD encounters divisions by potentially very small 
singular val ues, which is a strong source of potential numerical inaccuracies; but these can be 
eliminated jlOSll at low numerical cost. From a numerical point of view, tDMRG is not just a 
translation of TEBD, which came first, but an algorithm of its own, with strengths and weak- 
nesses. 

Both methods share the central feature that time evolution and truncation are intertwined: 
after each bond evolution, there is a truncation by SVD. By contrast, tMPS evolves all bonds 
first, and then truncates the entire state by compression of matrix dimensions SD — > D by SVD 
or iteratively. 

tMPS is the cleaner approach, but it can also be shown to be more precise. In fact, for 
real-time evolution it relates to tDMRG or TEBD exactly as iterative variational compression to 
compression by SVD, which implies that for small state changes (e.g. for very small time steps) 
the difference goes down, as the interdependence of truncations becomes less severe, there being 
only very benign truncations. That the above relationship exists can be seen from compressing a 
tMPS state not variationally, but by SVD only: 

Take to be right-canonical, and do a tDMRG/TEBD step on the first bond or tMPS steps 
on all odd bonds. The truncation is now to be carried out by SVD and my claim is that SVD does 
not see a difference between the two very different time-evolved states. On the first bond itself, 
all methods produce the same structure, but they differ on all other sites. Whereas 




(274) 




(275) 



cr 
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is 

l<^>tMPS = X 'i"'"''B"'T' . . . \CT}, (276) 
cr 

where the Z?-matrices come from the contraction with the time evolution bond operators. The 
SVDs on the first bonds are equivalent for both states provided both sets {B} and {B} generate 
sets of orthonormal states. This is indeed the case, because the B-matrices do this by definition, 
and the states generated by the S-matrices are related to the first set of orthonormal states by 
a unitary transformation (real-time evolution!). This observation of the equivalence of methods 
also holds for bonds further down the chain. 

Hence, the difference between the three algorithms becomes only visible at the level of vari- 
ational compression. 



7.4. How far can we go? 

In this section, I have described basic algorithms for the time evolution of pure and mixed 
states. There were two sources of error. One of them is the Trotter decomposition, which for an 
nth order decomposition generated an error (9(t"^^) for each time step t. As there are f/r time 
steps, the error will ultimately be 0{T"t), i.e. linear in time. This means it is only growing moder- 
ately in time and can be scaled down by smaller time steps and/or higher-order decompositions. 
This is common to all current methods |47. 48, 4^ S^, 51]. In fact, there are other methods of 
calculating matrix exponentials such as the Krylov method f 106] or lookahead procedures such 
as in 1107] , which reduce this error even more. In any case, it is not very worrisome in the long 
run. 

On the other hand, there is the error due to the truncation of the blown-up bond dimensions 
of the MPS after each time step. Th is err or is serious; early on it could be shown to lead to errors 
exponentially blowing up in time lll08ll . Yet truncation errors are only the symptom, not the 
fundamental problem: the real reason is that - following the Lieb-Robertson theorem - entan- 
glement S can grow up to linearly in time for an out-of-equilibrium evolution of a quantum state: 
S(t) < 5(0 ) + ct, where c is some constant related to the propagation speed of excitations in the 
lattice] 109]. This linear bound is actually reached for many quantum quenches, where a Hamil- 
tonian parameter is abruptly changed s uch t hat the global energy changes extensively. Both from 
D ~ 2^ and from a rigorous analysis jl lOtl it follows that in such cases the matrix dimensions 
will have to go up exponentially in time, D(t) ~ 2', or that for fixed matrix dimensions precision 
will deteriorate exponentially. 

Nevertheless, in many circumstances matrix size growth is slow enough that numerical re- 
sources are sufficient to observe the time-dependent phenomenon of interest: Time-dependent 
DMRG has been used extensively in the meantime and found to open completely new perspec- 
tives on the n on-equilibrium behaviour of strong ly correlated one-dimensional systems (to name 

a few: iscilslfmfiM fin flTi fill filing 



8. Time-dependent Simulation: Extending The Range 

Time evolution - whether it is done by TEBD, tDMRG or tMPS, to give the historical order - 
is fundamentally limited by the times that can be reached. The underlying reason is the (at worst) 
linear buildup of entanglement in time in an out-of-equilibrium quantum state, that translates it- 
self into an (at worst) exponential growth of bond dimensions D(t) if a given precision is desired. 
A "time wall" is hit exponentially fast. Can one push it further into the future? A similar issue 
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arises for finite-temperature calculations. While they are not necessarily about dynamics, seen as 
imaginary time evolutions they raise their own problems, regarding the T — > limit for static or 
thermodynamic calculations, and regarding dynamics, there in particular at high temperatures. 

A second issue that we have not covered so far concerns dissipative time evolution where 
we are not looking at a closed quantum system as in pure Hamiltonian dynamics but at an open 
quantum system. If the dynamics is Markovian (i.e. the bath has no memory, a highly non-trivial 
assumption in many cases), then the most general dynamics is given by the Lindblad equation. 
While it is easy to show that formally this can be simulated using MPS quite easily, in actual 
practice this is numerically involved and simpler schemes are highly desirable. 

In this section we will first consider attempts to extend the time range of simulation by dif- 
ferent schemes for evaluating the time evolution tensor networks. As we have already seen for 
simple examples like the evaluation of wave function overlaps, the order of contractions may 
hugely change the computational effort. In a second step, we will look at a prediction method 
that picks up on numerical raw data and extrapolates them very successfully over an order of 
magnitude, provided they meet a certain mathematical form, taking the case of finite temperature 
as an example. In a third step, we will look at an altogether different way of finite temperature 
simulations. In a last step, I will take up the issue of dissipative dynamics and show neat progress 
made in that field. 

8.1. Orders of contraction: light cones, Heisenberg picture and transverse folding 
Let us consider the calculation of the time-dependent expectation value 

imiOmt)} = <<Ale^"^'OPe-'^V>- (277) 

Starting with It//}, we evolve it up to time f, obtaining |i/'(f))- The expectation value then is 
calculated by sandwiching the two operators between (i/'(f)l and \if/{t)}, as discussed before. But 
we can represent this procedure also as the (approximate) contraction over a two-dimensional 
tensor network as shown in Fig.|50l which is then contracted line by line along the time direction, 
moving inwards. 

Assuming f - NAt and n MPOs per Trotter step (e.g. 2 in first order), we have a lattice of 
L X {2nN + 3) sites, i.e. of width L and odd height 2nN + 3. If we call r''-^' the tensor located 
on the site in row / and column j (like in a matrix), and if we label indices by up m, down d, left 
I and right r, and write r''-^' in analogy to MPOs with indices T^^'J^"''' ^ then we can identify, for 
example, for / - 1 (location of the bra state): 

ji\M\,d _ Amd* j[i,m,d _ Ai]d* j[\,m,d _ ALv* ,219.-) 

l,r l,r /,/■ /,r /,1 1,1 ^ 

where 1 < j < L. Similarly for / = 2nN + 3 (location of the ket state): 

j,[2nN+3,l]u,l _ .[Iju j,[2nN+i,j]ii,l _ .[j]u j,[2nA'+3,L]«,l _ (279) 

and on row / = nN + 2 (location of the operators): 

j,[nN+2,j]u,d _ 1 O"'"' on operator location / 
~ \ Su,d else 

In this row, horizontally the network is a product of scalars, hence the (1,1). On all other rows, the 
tensors T^''^^ are given by local MPOs such that j^'J^"''' - ^^'J"'"' on all rows nN+2 < i < 2nN+3 
(with the type a depending on the chosen decomposition) and T^''J^"''' = ^r[a]d,u* j.Q^g 
1 < ! < nN + 2, which correspond to the time evolution of the bra. 
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space 




space 

Figure 50: Two-dimensional tensor network contracted for time-dependent expectation values: the top and bottom lines 
represent (ifi\ and The two arrays of MPOs (indicated by brackets) represent e^'"' and e"''" in Trotterized form; I 
do not distinguish between different local MPOs such as identity operators which show up on some sites of the left- and 
rightmost columns. In the central line, we put identity operators (white squares) and the operators to be evaluated. The 
dashed line indicates the buildup of contractions in time direction. 
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8.1.1. Light cones in time evolution 

Considering the time evolution of bra and ket together in fact allows important simplifica- 
tions. In Fig.|5T|l have restored the alternating pattern of bond evolutions in a first order Trotter 
decomposition and explicitly marked the position of unit operators by white squares. We would 
like to calculate {0{t)), where the operator sits on site 2. Let us look at the last Trotter steps 
(rows 5 and 7). In row 5 there are several evolution operators e^'*"^' with corresponding operators 
g-i/iAr y g^j. jjjjg [jjgans that they cancel each other and can be replaced by unit operators, 

except in columns 2 and 3 because they have O interposed. If, in turn, we look now at rows 4 
and 8, there are now evolution operators cancelling each other in columns 5 through 8; the other 
ones do not cancel, as they sandwich non-identity operators. Like this, we work our way towards 
the bra and ket states, until no cancellation is possible anymore. The resulting tensor network 
shows a large degree of replacements of complicated tensors of bond (row) dimensions larger 
than 1 by identity tensors with bond dimension 1, which means that contractions become trivial 
and no compression is needed. There remains an algorithmic light cone of evolution operators 
that "see" the presence of a non-trivial operator O to be evaluated. Note that this algorithmic 
light cone is not to be confused with a physical light cone: if we send Af 0, the algorithmic 
light cone becomes infinitely wide for any fixed time f. Physically, the Lieb-Robinson theorem 
states that beyond a "light cone" of width x - 2ct, where c is some problem-specific "veloc- 
ity", correlations decay exponentially fast, as opposed to the hard cut imposed by a relativistic 
light cone. The physical light cone and the special decay of correlations is at the basis of very 
interesti ng algorith mic extensions of the MPS/tDMRG/TEBD algorithms of the last section by 
Hastingsl 105 . 1 16ll . which I will not pursue here. 

While this structure becomes more complicated if we look, e.g. at n-point correlators, we 
may look at a huge algorithmic saving, even though we have to pay the price that for different 
locations of operators, new networks have to be considered. 

What is the prize for calculating at different times, e.g. {0(t[)}, {0(t2)) and so on? This is a 
very natural question, as we might be interested in the time evolution of, say, some local density. 
If we do not use the light cone, then we simply calculate the contraction moving inwards, calcu- 
late some average, retrieve the stored result of the contraction up to the line with the operators, 
add more Trotter steps, contract, calculate some average, and so on. This is exactly what we 
have been doing all along. Of course, the light cone generated by the operator acting at time 
t2 > t\ is different from and larger than that generated by the operator acting at time ti . But if 
the Hamiltonian is time-independent, the larger light cone contains the smaller one at its tip. It 
therefore makes numerical sense to reverse the time evolution, and work from the future towards 
the past. But this corresponds to nothing else but a switch to the Heisenberg picture. 



8.1.2. Heisenberg picture 

Mathematically, the switch to the Heisenberg picture is nothing but a rebracketing: 

(o(/)) = m)\omt)) = <^ie+'^'oe-'^v) = {^mm, (28i) 

where we have introduced the time-dependent operator 

(9(f) = e^''"'de-'"'. (282) 

If we Trotterize the time evolutions present, we arrive exactly at the light cone structure of the 
last section, except that it has not been contracted yet with bra and ket; Fig.l52l 
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Figure 51: Light cone: Tlie diagrams are to be read from left to right, in the top and then the bottom line. On the top left, 
an operator is sandwiched between two time-evolved states; taking Trotter evolution into account, identity operators are 
marked by white squares. Working our way outward from the center towards the top and bottom, bond evolution operators 
cancel each other to identities, provided they do not sandwich the operator or (therefore) surviving bond operators. The 
result is a light cone of evolution operators, surrounded by numerically trivial identities. 




Figure 52: Time evolution of an operator O in the Heisenberg picture, translated to the MPS/MPO representation. Each 
symmetric set of layers around the operator corresponds to one time step (more precisely, part of a Trotter time step). 
The hght cone widens as a larger and larger section of the lattice is alfected by O. The identity operator (white square) is 
inserted for reasons of symmetry, but without explicit function. 
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This allows to set up time evolution in the Heisenberg picture 11521 l53h .Technicallv. one con- 
structs a spatially growing MPO from MPO-MPO-multiplications as encountered e.g. in the 
calculation of the action of fp-. If the current MPO of the time-evolved operator consists of local 
MPOs of the form <9^''[|a, bond dimension D (on sites / where it actually has to be consid- 
ered), and if the time evolution MPO (for e"''") reads W^'^^'^ with bond dimensions Dvy, then the 
operator reads after (part of) the time step 



{bi^\,ai-\,Ci-{),{bi,ai,Ci) 



bi-i.hi 



(283) 



with bond dimensions DD^. This operator is then compressed down to bond dimensions D as 
explained earlier for MPOs, essentially using the method for compressing MPS. 

This sets up a "conventional" time evolution: instead of a state, an operator in MPO form is 
subjected to time evolution by MPOs, and compressed to manageable matrix dimensions after 
each time step. We can basically recycle most of the algorithm. 

What are potential advantages and disadvantages of this formulation? First of all, the savings 
due to the algorithmic light cone are immediately incorporated. Second, we may hope that 
truncations are smaller: while the network contracted over is identical in the Heisenberg and 
Schrodinger picture, truncations in the Schrodinger picture do not take into account the operator 
and hence are less specific - one may enviseage that for "simple" operators like local density a 
lot of the fine structure of the state evolution is not really needed, and evolving the operator itself 
tells us which information is needed specifically for this operator 

A corresponding disadvantage is of course that calculations need to be redone for differ- 
ent operators, which in the Schrodinger picture may be evaluated whenever one likes, provided 
the time-evolved wave function is stored. Of course, here the corresponding advantage is that 
for different states one may evaluate whenever one likes, provided the time-evolved operator is 
stored. 

At the moment of writing it seems indeed that for simple operators the reachable time spans 
can be extended substantially, but I would find it hard to commit to some rule of thumb. 



8.1.3. Transverse contraction and folding 

Of course, the iterative build up of the state as it evolves in time appeals to our intuition about 
the world, but there is nothing that prevents us to contract the same network in the spatial direc- 
tion, i.e. column by column; as the order of contraction may influence efficiency quite strongly, 
maybe it helps. In order to recycle existing programs, one may simply rotate the current network 
by 90 degrees counterclockwise, and obtains a lattice of width {2nN + 3) and height L. If we 
continue to label tensors r'' '' by the vertical before the horizontal and in the row-column logic 
of a matrix, then tensors in the new lattice read 



l.r 



= T 



[L+l-j,i]r,l 
u.d ' 



(284) 



as can be seen from Fig.|53] Then we can contract again line by line. 

As it turns out, a simple rotation (or transverse contraction) does not extend the reachable 
timescale. It is by an additional /oWmg step that a strong extension of the timescale is possible 
ll54ll . The folding happens parallel to the new "time" (i.e. real space) axis, and halves the extent 
of the new "space" (i.e. real time) domain (see Fig.|54|i. Instead of sites 1 through 2nN -H 3 we 
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time [i] 



space [j] 
► 




rotation by 90° 
and reflection 
of „time'''-axis 



„space" [j] 
^ 



„time" [i] 



Figure 53: Rotating the lattice for code reusage: Assuming a space-time labeling with time / and space j (which 
after rotation refers to ficticious "time" and "space") tensor indices u,d and l,r exchange places as shown in the figure. 




folded „space" 



Figure 54: The rotated lattice is now folded on the "space" (former time) axis which has 2nN + 3 = L sites. Sites 
corresponding to same times come to cover each other (indicated by an ellipse); the line on which operators are evaluated 
at final time remains single at the bend. 



then have double sites 1 through nN + 2, where double site 1 comprises old sites 1 and 2nM + 3, 
double site 2 comprises old sites 2 and 2nN + 2; generally i comprises old sites / and 2nN + 4-i, 
up to the pair (nN + l,nN + 3). The site nN + 2 is special: as we are folding an odd number of 
sites, one site remains single. This is site nN + 2, which corresponds to the line that contained 
the operators to be evaluated at time t. On all the other sites, we fold tensors onto each other that 
correspond to "identical" timesteps, one forward and one backward in time. 

The expectation is that this folding of forward and backward timesteps leads to cancellations 
in entanglement buildup, such that larger times can be reached (the growth in D is not as fast). 

To simplify notation, we define L' = 2nN + 3 = 2^ + 1 . If we read the bottom end of the 
folding as an MPS, the folded state also starts with an MPS whose matrices are formed as 



for all 1 < / < ^ and 



^/K+US... ^ t^^'l^r = rl^f . (286) 



We have defined "fat" local states = |o",)|cr£,+i_,) on each site, except site f+l, where it re- 
mains of dimension d (for programming, one may of course introduce a dummy site). Similarly, 
we construct the new tensors for the folded MPOs. 

If we assume that the original lattice and the Hamiltonian acting on it were translationally 
invariant, at least for translations by an even number of lattice sites, we can write the contractions 
conveniently using a transfer operator. If we call the state on the first line (first "time" slice) of 
the folded lattice (corresponding to site L of the original lattice) and the one on the bottom 
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space 



Figure 55: Representation of the transfer operator E (dashed rectangle; note direction of action) in the unfolded, unrotated 
representation of the time evolution network. It repeats throughout the lattice by spatial translation by two sites, except 
on the sites with evaluation of an operator, where it is modified accordingly. 



line (last "time" slice) then (/ odd for simplicity) 

{>lrL\E^'-'^'^EoE(^-'-'^/^\ifr,) 

^""'■^ ^ <^.|/?<™l^.> ■ ^'''^ 

Here, we have introduced the transfer operators E and Eg on stripes of length { and width 2, 
as represented in Fig.|55](in unfolded, unrotated form for simplicity of representation). Eq is 
derived from E by inserting O instead of the identity at site £ + 1. 

This can be evaluated by iterative contractions and compressions for spatially finite lattices, 
but one can also take the thermodynamic limit. Let us assume an eigenvector decomposition of 
E as 

Note that E is not hermitian, hence |/) and (/| are not adjoint, but distinct right and left eigenvec- 
tors. From the biorthonormality of those, 

hm E^ = lim V Af\i}{i\ = A^\R}{L\, (289) 

/ 

where I have assumed that the largest eigenvalue /Iq is non-degenerate (which is usually the case) 
and changed notation to \R) and (L| for the associated right and left eigenvectors. 
We then obtain in the thermodynamic limit as expectation value 

. <^^, (290) 
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where Aq - {L\E\R). 2-point correlators would then be given by 



(diPj) = 



{L\EoE'-Ep\ 



]R} 



(291) 



where r is the number of transfer operators between the sites ; and j. 

In order to evaluate such expressions, we obviously need Aq, \R} and (L|. As the components 
of E are explicitly available, we can construct its transpose equally easily, hence reduce all to 
the determination of two right eigenvectors. As we are looking for the largest eigenvalue, the 
power method (iterative application of £ to a guess vector) will work. But one can equally read 
E as the MPO representation of some non-hermitian operator, and reuse iterative ground state 
search techniques, with two modifications: the search for the lowest eigenvalue is replaced by the 
highest eigenvalue and the conventional Lanczos algorithm has to be replaced by non-hermitian 
methods, either the biorthogonal Lanczos algorithm or the Arnoldi method. 

While the coding is more involved than for standard time evolution, the timescales reachable 
are extended substantially, factors 3 to 5 seem easily possible. 

8.2. Linear prediction and spectral functions 

Spectral functions S (k, cS) are among the most important theoretical and experimental quan- 
tities in r nany -body physics. While there are very accurate ways of calculating them directly at 
I12I1 . there is also an indirect approach, pioneered in llsol] . to calculate real-time 
real-space correlators like {St{t)S^.(Q)), and to carry out a double Fourier transform to momen- 
tum and frequency space. This approach has the advantage to extend to finite T seamlessly, but 
suffers from the limitations of reachable length and time scales. 

Of these, the limitations in time are much more serious, because of the rapid growth of entan- 
glement in time. The time scales reachable are mostly so limited that a naive Fourier transform 
gives strong aliasing or that one has to introduce a windowing of the raw data that smears out 
spectral information quite strongly. This limitation can be circumvented however at very low 
numerical cost by a linear prediction technique both at T = 0|117, 118] and T > O ll 1911 that 
extends reachable f and thereby greatly refines results in the frequency domain. 

For a time series of complex data xq, xi, . . . , x„, . . . , at equidistant points in time t„ - ntSt 

(and maximal time fobs := A^Af) obtained by DMRG one makes a prediction of xn+\,xn+2, 

For the data points beyond t - fobs, linear prediction makes the ansatz 



The (predicted) value x„ at time step n is assumed to be a linear combination of p previous 
values {jCn-i , . . . , x„-p]. Once the a, are determined from known data, they are used to calculate 
(an approximation of) all x„ with n > N. 

The coefficients a, are determined by minimizing the least square error in the predictions over 
a subinterval f„ e (fobs - ffit, fobs] of the known data (corresponding to a set Q = {n \ fobs - ffit < 
«Af < fobs)), i-C- we minimize in the simplest approach E = Y^nesi - x„\^. t^ = fobs/2 is often a 
robust choice to have little short-time influence and enough data points. Minimization of E with 
respect to a, yields the linear system 



p 




(292) 



i=I 



Ra = -r. 



(293) 
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where R and r are the autocorrelations 7? - 'ZjneQ.x*„_jX„-i and rj - Yj„eax*n-j^"- * I293I ) is 
solved by a = -R^^r. 

One may wonder why extrapolation towards infinite time is possible in this fashion. As 
demonstrated below, linear prediction generates a superposition of oscillating and exponentially 
decaying (or growing) terms, a type of time-dependence that emerges naturally in many-body 
physics: Green's functions of the typical form G(k,a>) - {co — et — 'E(k,a>)y^ are in time- 
momentum representation dominated by the poles; e.g. for a single simple pole at co - llii - irii 
with residue ci. Green's function will read G{k, t) - cie""^''"'''', and similarly it will be a super- 
position of such terms for more complicated pole structures. Often only few poles matter, and the 
ansatz of the linear prediction is well suited for the typical properties of the response quantities 
we are interested in. Where such an ansatz does not hold, the method is probably inadequate. 

To see the special form of time-series generated by the prediction, we introduce vectors 
x„ :- [x„, . . ., JCn-p+i]^ such that ( 12921 ) takes the form 



with 



A = 





Xyj + l 


— Ax^, 






-fl2 


—(23 


—a 


1 





■■■ 








1 


■■■ 










1 






(294) 
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with the fl, as the elements of the vector a found above. Prediction therefore corresponds to 
applying powers of A to the initial vector ic^. A (non-hermitian) eigenvector decomposition of A 
with eigenvalues a, leads to 



^[A^xlv]!^^^'^' 



(296) 



where coefficients c, are determined from x^ and the eigenvectors of A. The eigenvalues a, 
encode the physical resonance frequencies and dampings. The connection is given as ff, 

piai,A/-);,A/ 



. Spurious \ai\ > 1 may appear, but can be dealt withf 1 19]. 
At r = 0, critical one-dimensional systems exhibit power-law decays in their time-dependent 
corre lators. The superposition of exponential decays is then taken to mimic these power-laws 
1 1 ITIl . At finite temperatures, time-dependent correlators S(k,t) decay typically exponentially 
for large times (due to thermal broadening), making linear prediction especially well-suited for 
this situation. This is also close to typical experimental situations, like inelastic neutron scattering 
off one-dimensional magnetic chains. 

As example, let us consider a field-free Heisenberg antiferromagnet with P - (XF-chain) 
and J~ - \. The former case allows for an exact analytical solution. It turns out that prediction 
allows to extend time series S {k, t) by over an order of magnitude without appreciable loss of 
precision. In frequency space, this corresponds to extremely high-precision spectral lineshapes 
(Figurel56b. 

As the dispersion relation of the XK-chain is just a simple magnon line, its self-energy 
structure is very simple, hence the prediction method easily applicable. As a more demand- 
ing example, we consider the spinon continuum of an isotropic 5 = 1/2 chain; Fig.|57l In the 
zero-temperature limit, results agree extremely well with Bethe-ansatz results (where remain- 
ing differences are hard to attribute: the Bethe ansatz here can only be evaluated approximately 
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frequency to 

Figure 56: Lines and dots represent exact analytical and numerical solutions for the lineshape of the spectral function 
S*^(k, cj) of an Xy-chain at temperatures fS = 10 and yS = 50 (broad and narrow Uneshapes) for (from the left) k = rr/S, 
k = nj'i, k = 'in I A. The dashed lines are the shapes that would optimally be extracted from the /3 = 10 simulation without 
prediction using some windowing of the raw data before Fourier transformation. Adapted from Ill9ll . 

lll20ll '). At finite temperatures, simulations at different precision indicate that results are fully 
converged and essentially exact. This lets us expect that this method will be a powerful tool in 
e.g. simulating the results of neutron scattering experiments. 

8.3. Minimally entangled typical thermal states 

Simulating thermal density operators for the calculation of static and dynamic properties 
works very well. Theoretically, there are no limits to this method. In practice, one encounters 
several limitations. On the one hand, simulations become difficult at very low temperatures 
r ^ 0. In this limit, the mixed state living on the physical system P will evolve towards the pure 
state projector on the ground state, = \4'<x>)p p{4'ai\ (here I use |i/'oo)p as the ground state of the 
physical Hamiltonian H, to refer to y6 = oo as in the purification section). But in this limit P is not 
entangled with the auxiliary system Q anymore, and we simulate a product of two pure states: 
assume for simplicity that the ground state energy is set to 0. Consider now the reduced density 
operator for the auxiliary system. Up to an irrelevant norm, 

= lim Trpe-^«/Vo)(<Ao|e-^^^' = p<-Aoo|iAo><-AoI«Ac«)/>, (297) 

because in the limit — > oo the trace reduces to the ground state contribution, Tr^ — > /)(i/'oo| ■ 
|iAco)/)- With the y6 = purification of the density operator, |i/'o) = d^^^^ Ynj- \(t)p\(t)q, the last 
expression reduces, again up to an irrelevant norm, to 

= Z 'A^'l'^>e e<o-'l<A^' = \^^)q e<'Aoo|, (298) 

crcr' 

where the i}/^ are the expansion coefficients of the ground state. Hence, the zero temperature 
purification is a product state 

= l^oo)p|iAoo>e, (299) 

where the latter state is just the physical ground state defined on the auxiliary state space. As- 
suming that it can be described with sufficient precision using matrix dimension D, the product 
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Figure 57: Spectral function of tlie isotropic Heisenberg cliain at two mome nta, and four different temperatures. At 
T = 0, Betlie ansatz (B.A.) and numerics agree extremely well. Adapted from Ill9ll . 



will be described by matrices of dimension D^. Effectively this means that our algorithm scales 
with the sixth instead of the third power of the characteristic matrix dimension for the problem 
under study. On the other hand, and this is the issue prediction has tried to address, we en- 
counter long-time simulation problems in particular at high temperatures T — > oo: many states 
contribute at similar, but not identical, weight and MPS are not efficient at encoding this wealth 

of contributing state s. 

As White | lOOl 121 1 has pointed out, one can avoid the purification approach entirely by 



sampling over a cleverly chosen set of thermal states, the so-called minimally entangled typical 
thermal states (METTS). This approach has already been shown to alleviate strongly the first 
limitation, while not much is known yet about the second limitation. 
A thermal average is given by 

(A) = ixre-^^A = i ^ e-^^"<«|A|«), (300) 

n 

where 1 have chosen, like all textbooks do, the energy representation of the thermal density op- 
erator. As already pointed out by Schrodinger many decades ago, this is mathematically correct, 
but unphysical in the sense that real systems at finite temperature will usually not be in a statis- 
tical mixtures of eigenstates, as eigenstates are highly fragile under coupling to an environment. 
But the choice of the basis for taking the trace is arbitrary, and one may also write 



(A) 



J ^ </ie-^«/2^e-^^^2|/) ^\Yj mmmm) ood 



where {|/)) is an arbitrary orthonormal basis and |0(O) - P{iY^I'^er'^"'^\i). With P{i) - </|e"^^|/), 
we recognize to be normalized. It is easy to see that Y^iPii) - -Z, hence the P(i)/Z are 
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probabilities. One can therefore statistically estimate (A) by sampling \(p{i)) with probabilities 
P{i)IZ and average over ((^(i)|A|<^(/)). 

Several questions arise before this can be turned into a practical algorithm. How can we 
sample correctly given that we do not know the complicated probability distribution? Can we 
choose a set of states such that averages converge most rapidly? Given that an imaginary time 
evolution will be part of the algorithm, can we find a low-entanglement basis {\i)} such that time 
evolution can be done with modest D, i.e. runs fast? 

To address these issues, White chooses as orthonormal basis the computational basis formed 
from product states, 

= 1/1)1/2)1/3)... (302) 

classical (unentangled) product states (CPS) that can be represented exactly with D - \. The 
corresponding states 

l^^(O) = -4=^'^"'^\i) (303) 

are so-called minimally entangled typical thermal states (METTS): The hope is that while the 
imaginary time evolution introduces entanglement due to the action of the Hamiltonian, it is 
a reasonable expectation that the final entanglement will be lower than for similar evolutions 
of already entangled states. While this is not totally true in a strict mathematical sense, in a 
practical sense it seems to be! Compared to purification, this will be a much faster computation, 
in particular as the factorization issue of purification will not appear. 

In order to sample the !</>(/)) with the correct probability distribution, which we cannot cal- 
culate, one uses the same trick as in Monte Carlo and generates a Markov chain of states, 
I'l) ~* \h) — > l'3) — > ... such that the correct probability distribution is reproduced. From 
this distribution, we can generate \^{ii)), |<^(/2)), \<p{i-i)), . . . for calculating the average. 

The algorithm runs as follows: we start with a random CPS |/). From this, we repeat the 
following three steps until the statistics of the result is good enough: 

• Calculate er^^^^\i) by imaginary time evolution and normalize the state (the squared norm 
is f (/), but we won't need it in the algorithm). 

• Evaluate desired quantities as (0(/)|A|0(/)) for averaging. 

• Collapse the state |0(/)) to a new CPS |/') by quantum measurements with probability 
p{i — > /') = K/'|0(/))p, and restart with this new state. 

Let us convince ourselves that this gives the correct sampling, following floS). As the |0(/)) 
follow the same distribution as the |/), one only has to show that the latter are sampled correctly. 
Asking with which probability one collapes into some |j) provided the previous CPS |/) was 
chosen with the right probability P(/)/Z, one finds 

^ ^p(/ - ;■) = 2 ^\<jmt = 2 ^^m-'^'I'm' = ^01e-^^|;) = (304) 
(' / / 

This shows that the desired distribution is a fixpoint of the update procedure. It is therefore 
valid, but it is of course sensible to discard, as in Monte Carlo, a number of early data points, to 
eliminate the bias due to the initial CPS. It turns out that - after discarding the first few METTS, 
to eliminate effects of the initial choice - averaging quantities over only a hundred or so allows 
to calculate local static quantities (magnetizations, bond energies) with high accuracy. 
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While we already know how to do an imaginary time evolution, we still have to discuss the 
collapse procedure. As it turns out, the structure of MPS can be exploited to make this part of 
the algorithm extremely fast compared to the imaginary time evolution. 

For each site /, we choose an arbitrary li-dimensional orthonormal basis {|(t,)), to be distin- 
guished from the computational basis {|cr,)). From this we can form projectors P""' - |a",)((T,| 
with the standard quantum mechanical probability of a local collapse into state |a",) given by 
Pa-i - (i/'l^'^'I'A)- If we collapse li/r) into the CPS |(/'') - |<5'i)|5"2) . . . Iitl), the probability is given 
by {ij/\P°'' . . .P°'''\if/} = as demanded by the algorithm. After a single-site collapse, the 

wave function reads 

m ^ p-^'/^P'^liP), (305) 

where the prefactor ensures proper normalization of the collapsed state as in elementary quan- 
tum mechanics. To give an example, for S = ^ spins measured along an arbitrary axis n, the 
projectors would read 

pT,„i„ = i + n ■ S,-. (306) 

Such a sequen ce of local measurements and collapses on all sites can be done very efficiently, 
as pointed out by 1 100], if one exploits two features of MPS and CPS: (i) local expectation values 
can be evaluated very efficiently if they are on explicit sites (in DMRG language) or on sites 
between left- and right-normalized sites of a mixed-canonical state (in MPS language) and (ii) 
after the collapse, the resulting state is a product state of local states on all collapsed sites and 
ffie uncollapsed remainder. 

Assume that is right-canonical (with the relaxation that the normalization on site 1 is 
irrelevant). Then the evaluation of pg-^ - {ip\P°'' \^) trivializes because the contraction of the 
expectation value network over sites 2 through L just yields 5^1, a' ■ Hence 



(307) 



This expression looks very specific to the first site (because of the open boundary), but as we will 
see it is not! 

Once the probabilites for the collapse on site 1 are calculated, one particular collapse is 
chosen randomly according to the distribution just generated, \&\). The state after collapse will 
be of the form \&i )|iA^rest), hence a product state. Therefore, the new matrices (which we call A'^' ) 
on site 1 must all be scalars, i.e. D - \ matrices. From \&) - 2^ \cr){cr\&) - 2^ \o')A°' they are 
given by 

Al\^{(T,\(ri)- (308) 

it is easy to see that left-normalization is trivially ensured, hence the labelling by A. But this 
change in the dimension of A'^' means that B"'- has to be changed too, namely 

(T,,ai 

As the label ixi takes a definite value, it is just a dummy index, and the row dimension of M°"- is 
just 1, like for ffie matrices on ffie first site. Hence, Eq. (13071 ) generalizes to all sites, and the most 
costly step is the update of B°'', which scales as D'^d^, but not as D^, as time evolution does. 
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To see the substitution, we express P°'^ as an MPO, W°''°'i - {(Ti|cri){iTi|crj). Hence, the 
collapsed \tp) reads 



cr a-' 



.02 



(B'^-V..B'^%,i|o->, 



which yields the substitution. 

A few more comments are in order At each site, the measurement basis can be chosen ran- 
domly, and in order to obtain short autocorrelation "times" of the Markov chain, i.e. high quality 
of the sampling, this is certainly excellent, but also much more costly than collapsing always into 
the same basis, which however generates ergodicity problems. The proposal is to switch alternat- 
ingly between two bases where for each basis projectors are maximally mixed in the other basis 
(e.g. if we measure spins alternatingly along the x- and z- (ory-)axis). Autocorrelation times then 
may go down to 5 steps or so jlOOll . For estimating the statistical error, in the simplest cases it 
is enough to calculate averages over bins larger than the autocorrelation time, and to look at the 
statistical distribution of these bin averages to get an error bar. 

Intriguing questions remain, concerning both the potential and the foundation of the algo- 
rithm: how well will it perform for longer-ranged correlators, as needed for structure func- 
tions? Dynamical quantities can be accessed easily, as the time-evolution of the weakly en- 
tangled METTS is not costly - but will the efficiency of averaging over only a few "typical" 
states continue to hold? 



8.4. Dissipative dynamics: quantum jumps 

Dissipative (i.e. non-Hamiltonian) dynamics occurs when our physical system A is coupling 
to some environment B such that A is an open quantum system. This is a very rich field of 
physics, so let me review a few core results useful here. The time evolution of the density 
operator of the system can always be written in the Kraus representation as 

PAit) = ^i(OPA(0)£f (f), (310) 

j 

where the Kraus operators meet the condition 

2£f(f)£i(0 = /A. (311) 

j 

If the dynamics is without memory (Markovian), it depends only on the density operator at an 
infinitesimally earlier time, and a master equation, the Lindblad equation, can be derived. In 
the limit df — > 0, the environment remains unchanged with probability pQ — > 1 and changes 
(quantum jumps) with a probability linear in df. If we associate Kraus operator Ef^it) with the 
absence of change, a meaningful ansatz scaling out time is 

£^(df) = /a + 0(df) E{(dt) ^ yfdtL{ (;>0) (312) 



or more precisely 



£°(dO = /a + {Ka - ii/A)d?, 
100 



(313) 



with two Hermitian operators Ka and Ha ■ The normaUzation condition of the Kraus operators 
entails 

»o 

These ansatzes allow to derive a differential equation from the Kraus evolution formula, which 
is the Lindblad equation 



^ = -i[H,p] + g {LjpLj-^- - ^{L^U^plj , (315) 

where I have dropped the indices A. Indeed, in the absence of quantum jumps (j = only), 
one recovers the von Neumann equation. At the price of non-hermiticity, this equation can be 
simphfied. If we introduce i/eff H + iK, then the last term disappears and we have 

^ = -i[i?eifP - Mrff] + 2 ^P^'- (316) 

The simulation of Lindblad equations is possible quite easily in the MP S form alism, in par- 



ticular using MPOs |51], but also in the form of a superoperator formalism 1112411 . The problem 
with this approach is that it is numerically more costly compared to the Hamiltonian evolution 
of a state. A very attractive alternative, which allows maximal reusage of available pure state 
codes, has been proposed by 1 125], which combines pure state time evolution with the method 
of quantum trajectories. 

The method of quantum trajectories has been widely applied in quantum optics f 126']. Instead 
of using the Lindblad equation directly (which takes into account both the probabilistic distri- 
bution of initial states through p(0) and all possible sequences of quantum jumps), the quantum 
trajectory approach samples over the distribution of initial states, and for each of this sample 
states carries out a pure state time evolution where random quantum jumps occur at random 
times. They are chosen such that if one averages physical quantities over this distribution of 
time-evolving states, the result of the Lindblad equation is recovered. Let us ignore the sampling 
over initial states, assume that it is always the same, and instead focus on the conceptually more 
difficult averaging over quantum jumps. 

The algorithm then proceeds by generating N quantum trajectories in a time interval [0, T] 
(where T is the final time of the simulation) as follows: 

• Generate a starting state |i/'(0)); it either samples the t - density operator correctly or is 
simply always the same, depending on the physical problem. 

• Choose a uniformly distributed random number po in [0, 1]. 

• Carry out, using one of our pure state time evolution methods, the time evolution of |i/'(0)) 
under i/gff ■ As the effective Hamiltonian is non-hermitian, the norm of the state will de- 
crease over time. Stop the time evolution at time fi, which is defined by {i//{ti)\il/{ti)) - po; 
this is the time of the first quantum jump. Note that if T < fi, our simulation stops at T 
and we have a trajectory without jump, and we normalize the final state. 

• To carry out the quantumjump at t[, we calculate 

pj^miW^Ljmi)} pj^-fi- (j>o) (317) 

2j >>o Pj 
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and choose a / according to the normaUzed probability distribution {pj}. 

• We carry out this jump and normalize the state, 

IW> = 4^^. (318) 

• After this, we continue with finding a new po, from which time evolution of |i/'(f|)) with 
^eff generates f2, the location of the second quantum jump, and so on, until T is exceeded. 

Physical quantities up to time T are now averaged over the TV quantum trajectories that have 
been generated. The correct probabilities are produced if all states are normalized at all times; as 
this is not the case in the algorithm, norms at say time t have to be taken into account. Obviously, 
a careful analysis of convergence in TV — > oo has to be carried out, but it seems that fo r a s mall 
number of jump operators, even a few 100 trajectories may give highly reliable results f 125']. 

The observation that this sampling reproduces the dynamics of the Lindblad equation is part 
of the standard literature on quantum trajactories. The proof can be done in two steps, which I 
just sketch here. In a first step, one considers fixed time steps dt, and calculates probabilities for 
nojump vs. jump y in this time interval (p^ - At{i^{t)\V'^ V\ij/{t)), p^) - l-TjjPj)- One then either 
time-evolves under i/eff over df and normalizes, or does the jump and normalizes, according to 
the generated distribution. One can show that this reproduces the Lindblad equation. In a second 
step, one shows that the distributions of quantum jumps generated in this way and the one we 
use in the algorithm are identical. 



9. DMRG and NRG 

9.1. Wilson 's numerical renonnalization group (NRG) and MPS 

Wilson's Numerical Renormalization Group (NRG) 63, ill] originates in attempts to 



explain why metals with a small concentration of magnetic impurities exhibit a non-monotonic 
behaviour of resistivity. It was found that an adequate minimal model is provided by 

Ha^Y, ''^l^k. + Z ^^(f^^ka + h-c.) + U(h„ - l/2)(hf^ - 1/2). (319) 

kcr kcr 

This single-impurity Anderson model contains an impurity site that can be occupied by up to 
two electrons (operators f^) with on-site repulsion U and which couples to a conduction band 
(operators c|,^) with energy dispersion through some hybridization function Vk- 

In order to make it tractable, one changes from momentum to energy representation, as- 
suming that only low-energy isotropic i-wave scattering matters, and introduces logarithmic 
discretization: the band is represented by band segments of an energy width that decreases ex- 
ponentially close to the Fermi energy €f. This accounts for the observation that the decisive 
feature of quantum impurity physics, namely the appearance of a very narrow resonance peak at 
the Fermi energy in the local impurity spectral function, is linked exponentially strongly to the 
states close to the Fermi energy. Logarithmic discretization is however also required to make 
NRG work at all on a technical level! 
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After further manipulations, for which I refer to lf59l.l60ll . the Anderson Hamiltonian is finally 
mapped to a semi-infinite chain of non-interacting sites with the exception of the first one: 



H = f/(n/t - l/2)(n/x - 1/2) + f_i ^(/J^oc. + h.c.) + ^ f„((JI^4+i,^ + h.c), (320) 

a- (r,n=0 

where the Jo- are fermionic operators. The crucial point is that the t„ decay exponentially, t„ ~ 
A"", where A is the shrinkage factor of the energy bands in the logarithmic discretization, usually 
a value of the order 1.5 to 2. This is obviously a model that is amenable to our methods, e.g. a 
ground state search - as the hoppings decay exponentially, we will not have to consider a truly 
infinite chain. 

NRG builds on the observation that the exponential decay leads to a separation of energy 
scales: assuming we know the spectrum of the partial chain up to some length, all remaining 
sites will only make exponentially small corrections to it because of the exponentially small en- 
ergy scales further down the chain. Finding the ground state (and more generally the low lying 
spectrum) is now achieved by iterative exact diagonalization: assume that we have an effective 
D-dimensional eigenspace for some left-end part of the chain. Then the next-larger chain has 
state space dimension dD - 4D; in order to avoid exponential growth, we have to truncate down 
to D states. The NRG prescription is to diagonalize that system and to retain the D lowest-lying 
eigenstates. Starting out from very short chains that can still be done exactly, this procedure 
resolves the lowest-lying states exponentially well and is justified by the separation of energy 
scales: the decision which states to retain at some step would not be drastically changed with 
hindsight, as all further sites in the chain interact at much smaller energies. The obtained eigen- 
spectra at different energy scales (chain lengths) can then be used to extract RG flow information 
or calculate thermodynamic or dynamic quantities for the impurity problem. 

Given that the building block A°" of an MPS can be interpreted as encoding a decimation step 
upon growing a block by one site, irrespective of the decimation prescription, it is immediately 
obvious that NRG, like DMRG, can be seen as operating on MPS i62tl . This closes a historical 
loop as in fact the analysis of failures of NRG naively applied to Heisenberg and Hubbard models 
gave rise to the development of DMRG. A NRG state would look like 

|flf)= 2 (A'^^ ...A''%\cTi...at}. (321) 

(Tl,...,crf 

At each length {, we get a spectrum of D states. 

Given that DMRG is variational over the MPS ansatz space, it is reasonable to expect that at 
least some improvement must be possible over the NRG method. In fact this is the case |62]; in 
the next section, I am going to discuss some improvements which are already firmly established 
and others which are more speculative, i.e. where benchmarking on relevant complex problems 
is still lacking. 



9.2. Going beyond the numerical renormalization group 

In fact, considering an MPS formulation of NRG helps even without resorting to the connec- 
tion t o var i ation al methods like DMRG, as exemplified by the strict enforcement of certain sum 
rules 0127 . 128 1. but this is outside the topic of this review paper. 

What we can do more, however, is to subject the final MPS construction generated by NRG 
to DMRG-hke sweeping. This will somewhat improve the quality of the ground state, but above 
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Figure 58: Impurity spectral function for a Kondo model in an external field. Thin veilical lines indicate the energy inter- 
vals: below a certain energy, logarithmic discretization is turned linear. Variational MPS calculations (with and without 
deconvolutio n) re veal a peak missed by NRG, where the peak position is in very good agreement with a calculation by 
Rosch et al. Il29ll combined with a perturbatively found peak shift. Taken from (631 ■ 



all, the truncation procedure for high energies (short chains) will learn about truncation at low 
energies and vice versa. As opposed to NRG, there is now a feedback between energy scales. In 
that sense, NRG for an impurity problem is a similar conceptual step as the warm-up procedure 
infinite-system DMRG provides for variational finite-system DMRG. 

For logarithmic discretization, energy scale separation is big enough that this effect is minor 
and for a simple single impurity problem with a focus on the Abrikosov-Kondo-Suhl resonance 
the ultimate improvement is very limited, as NRG is geared to describe this feature optimally. 
The essential point is that energy scale separation can now be abandoned altogether due to feed- 
back, hence also logarithmic discretization, and we may choose a more fine-grained resolution 
of the energy band wherever it is physically suitable. This could find a variety of applications. 

In one application, variational calculus over MPS was applied to an impurity problem in an 
external field. The external field leads to a splitting of the peak into two spin-dependent ones, 
shifted above and below the Fermi energy. In Figure |58] we consider one of these peaks, using 
three techniques, NRG, an analytical approachl 129, 130], and variational MPS (DMRG) calcu- 
lus. NRG due to logarithmic discretization focuses on ep and does not see the field-dependent 
peak at all. Relaxing logarithmic discretization and providing sufficiently fine energy intervals 
around the expected peak positions away from ep the shifted resonance can be resolved clearly 
and even in very good agreement with analytics. 

A second interesting application of this could be to replace NRG as an impurity solver in 
the context of the dynamical mean-field theory (DMFT) IJJl, 132, 133, 134, IJJJ. In that case, 
information beyond the metallic resonance at the Fermi energy is required such that improving 
spectral resolution on other energy scales would be highly desirable. 

As the semi-infinite chain is non-interacting but on the first site, one can think about un- 
folding it into an infinite chain of spin- 1/2, with the impurity at the center and the presence or 
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absence of spin-up or spin-down fermions corresponding to the 2 spin states, the left half of 
the c hain corresponding to the spin-up fermions and the right half to the spin-down fermions 
1 13411 . Similar energies are now no longer grouped together, but in a DMRG-like approach this 
does not matter anymore! The intuition that spins that interact only through the central impurity 
might be essentially unentangled is corroborated by actual calculations. This is important as this 
means we will not pay a strong price by increased matrix dimensions. On the contrary: if in the 
NRG approach we are essentially looking at two uncoupled spin chains parallel to each other, 
this means that the corresponding MPS has dimension 0{D^) if the spin chain has dimension D. 
We can therefore expect that a NRG calculation with state number D can be replaced by a faster 
DMRG calculation with a state number 0{ VD). 

Beyond this speedup, unfolding can of course also b e app lied if the impurity couples to multi- 
ple bands, where NRG becomes exponentially complex ll62ll . The central site, of course, remains 
the same, and its numerical treatment can become extremely costly, such that new strategies have 
to be designed for that site. M uch work r emains to be done here, but first interesting follow-ups 
on these ideas have been made 1 122, 1231. 



10. Infinite-size algoritlims 

70.7. Reformulation of infinite-system DMRG in MPS language 

After the extensive discussion of finite-system algorithms, let us now reformulate infinite- 
system DMRG entirely in MPS language. It is convenient to label local states a bit differently 
to account for the iterative insertion of sites; we call the states |cr^)|cr^) . . . |cr^)|cr^) . . . Ictj )|cr^). 
Moreover, it will be very useful to give two labels to the matrices A and B, because the link 
between the matrices and the site on which they were generated will disappear 

Starting from blocks of size (i.e. a chain of 2 sites), the ground state wavefunction is - 
2^o-s *P°^°"i'|cr^)|cr^). Reading '^'^'1"'" as matrix ^.b, it is singular-value decomposed as 

^F' = t/iAlHyJ . From this we read off 

a'I!:: - (U,)^,, <!f = (Vj),,^. (322) 

A and B inherit left and right-normalization properties from U and V^, and the state takes the 
form 

l^i) = 2 Al'l-^Al'lBl'l-flcrffrf). (323) 

If we now insert two sites, and minimize the energy with respect to H, we obtain 

|^2>= 2 A[il'^>F^'^?B[')'^?|crfc4(rfcrf), (324) 

where each ^'°i°'2 is a matrix with dimensions to match those of A and B, implicit matrix multi- 
plications A'i'B assumed. Reshaping this set of 'I'-matrices into one, 

(*i'2)(afa^),(afo-f) = (^'^°"- ,af , (325) 

SVD gives = f/aA™ from which we can form 

105 



i i i i»i i ; i 

Figure 59: MPS structure generated by infinite-system DMRG at step C: a string of left-normalized A, a string of right- 
normalized B, joined by a diagonal singular value matiix A''^' . Note that structurally the central unit does not repeat. 

such that 

1^2)= 2 Al'l-^APHAraBPHBliKl^^^B^B). (327) 

At the ^th step, the wavefunction will read 

1^^,)= ^ A[iJ'^.\..AM'^'AMBM<...BnK|(^...(^(^fi...(ri*) (328) 

a^...a^a-f...a-^ 

and look like in Fig.l59l 

Of course, at each step we discard the smallest singular values and their associated singular 
vectors once matrix dimensions exceed D, which is nothing but the density-matrix based trunca- 
tion in the original formulation. At each step (new chain length) we can write down H for that 
length as an MPO and do the energy minimization. Other operators find similar representations 
as in the finite-size case. 

Let me briefly go through the reformulation of this algorithm in the FA-notation. In the 
first step we simply rename A and B into F, in line with the translation of boundary sites in the 
finite-system case: 

1^1) = 2 FliKAl'lrlil-^l^rftrf). (329) 

We then minimize in 

1^2)= 2 Fl'l'^'P'^'^rl'l-^lfrfcr^trfcrf), (330) 

and decompose it - as before - into A'^^°^ A'^^B'^^"^. Now we define (and due to the labelling, 
there is a slight change for the B-matrices compared to the finite-system setup) 

a['1f1=a'^^^' r<.?.[i]^5i2H (331) 

ab ah ab h ah ^ 

to arrive at 

1^2)= 2 F['l'^'A['lFPHAPlFPl-2A[ilF['l-?|crf(r^crfcrf), (332) 

as represented in Fig.l60l 

We can now ask two questions: (i) in DMRG, finding the ground state by an iterative solver 
like Lanczos is the most time-consuming part. Can we find a speed-up by providing a good initial 
guess? In finite-system DMRG the MPS formulation automatically yielded White's prediction 
method, whereas attempts at speeding up infinite-system DMRG have been made in the past, 
meeting with mixed success [1 36. .137il . (ii) Can we use the information at the chain center to 
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Figure 60: MPS structure generated by infinite-system DMRG at step ( in tlie FA-notation: F and A matrices alternate, 
and a possible identification of a (repetitive) two-site building block is given. 



build a translationally invariant state (up to period 2) in the thermodynamic limit, find its ground 
state or evolve it in time? 

The answer is yes to both questions, and builds on the identification of a two-site repetititve 
structure in the states. As opposed to the A, B-notation, where the central unit does not repeat 
itself even in the thermodynamic limit, it is very easy to read off a two-site repeat unit in the 
FA-notation, given by 

A[f-i]r[fl-?A[flr[^l-'. (333) 
Using the translation rules it can be translated into the A, Z?-language: 

AM'^AmBW-f(A[f-il)-i. (334) 

This result can also be obtained directly from the A, B notation, but the argument is more compli- 
cated than in the FA notation. It is of course to be understood that repeating these state fragments 
does not generate the state they were taken from; the claim is just that in the thermodynamic limit 
^ — > oo, when all sites are created equal, this repetition can come close. In any case, they are an 
educated guess about what the state will look like! 

We will now put this state fragment to multiple use, first on finite systems generated by 
infinite-system DMRG and then on thermodynamic limit states, both in the context of ground 
state searches and time evolutions. In the former case, it will provide a highly eflicient guess 
for the next quantum state; the evaluation of observables on this state proceed exactly as in the 
other finite systems. In the second case, both ground state and time evolution algorithms can be 
formulated (iDMRG and iTEBD), which however necessitate both an (identical) analysis of the 
issue of orthonormality of states in the thermodynamic limit. 



70.2. State prediction in infinite-system DMRG 

The identification of the "unit cell" of the state allows us to define a good initial guess for 
infinite-system DMRG l,65i1 , which avoids all the problems encountered by previous authors and 
leads to a dramatic speed-up even for small chains, where the underlying assumption that the 
chain center is representative of the physics of the thermodynamic limit state is certainly wrong: 
in order to grow the chain, we simply insert one unit cell, even though for small chains the idea 
that the state is just a repetition of these unit cells is not well verified - but even then so much 
better than a random guess. Starting from 

1^,) = ^Al'J"^ . . .a[^-iJ<.(AM-?A[^1bM-?[A[^-"]-1)A[^-i1b[^-'1-?-. . . .fil^'^V), (335) 
cr 

where the repeat unit has been bracketed out, the guess will then read 
cr 
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(A[^KAMBm<i[A[f-ii]-i)(^m<,/^K]gm^rf|-^[f-i]pi-)^ 

A[^-"B[^-"<....Bl'l-V) (336) 



or, multiplying out, 

l^s-s^ ^ ^Al'J^? . . .aM'^AMbM'^..[A[^-"]-'aM<.AMbM'^B[^-'I-?-. . . (337) 
cr 

In this ansatz, we can now identify a guess 

"P^^r"' = Al^JfiM^. [A1^-'1]-'aM<. AM. (338) 

From this ansatz, we can then iteratively find the that minimizes the energy in the 

infinite-system DMRG framework, generating from it At'^^'^'^+i , A''^^'', and Z?l^+'l""f+i . 

Alternatively, the ansatz can be brought into in a more elegant form. At the moment, B- 
matrices show up on the A-side of the lattice and vice versa. But we can exploit our ability to 
canonize MPS, and canonize Al^lfit^Ki by SVD to Al^+'Ki Ajfl, where A is derived from U 
and A from DY' in the way described before (A^B^^ = Maa,b = Uao-,kDk{y"^)k,b = Kk^kh)- 
Similarly, we do a canonization from the right on A'^^°^+i A'^'^' to obtain aJ^^B'^^^'^'^+i, where B is 
from y ' . Then we have an ansatz 

l^guess^ = ^Al'l'^' . . .Al^+U'^f-'Alt'jB"^""^*' • . .B'^^'lo-), (339) 
cr 

where 

A^t!]=AM[A[^->l]-'AM. (340) 

From this ansatz, we can then iteratively find the A'^^'^ that minimizes the energy, slightly mod- 
ifying the minimization part of variational MPS for a single site. In general, the result will not 
have the diagonal form resulting from an SVD, because Ar and Ki are not diagonal to begin 
with. But an SVD on it yields two unitaries that can be absorbed into the neighbouring A and 
B without affecting their normalization properties, such that the final A''^^'' is diagonal. In this 
form, the algorithm can be represented as in Fig. 1611 

As shown by McCulloch i65ll . this prediction leads to a dramatic speedup of infinite-system 
DMRG which complements nicely prediction algorithms of finite-system DMRG: the overlap 
between the predicted and calculated state often approaches unity up to 10"'" or so! 

10.3. iDMRG: variational ground state search in the thermodynamic limit 

Using the ideas of the preceding sections, it is very simple now to turn infinite-system DMRG 
into a performing and precise algorithm, called iDMRG, referring to the thermodynamic limit 
version of DMRG: 

• Set up an infinite-system DMRG algorithm. 

• Add the prediction procedure of the last section to the minimization algorithm. 

• Run the modified algorithm until convergence is achieved. Convergence of the wave func- 
tion to the thermodynamic limit can be judged by considering the relationship Eq. (1121b 

^A[^l>[^U[^J-^'=pf 1, (341) 
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Figure 61 : Summaiizing representation of tlie infinite system algorithm including prediction. At each step, two more sites 
are introduced, with an ansatz (white) calculated from the last iteration, leading to a new singular value decomposition 
(black). 



where p^^^ - a'^^A'^^' and p^'^ = ^['"-i]^[''-iJt ^j-g j-jjg reduced density operators of the 
left half of the system. Note that while this relationship holds always between reduced den- 
sity operators in the same finite system, here they originate from systems of two different 
lengths 2(i - 1) and 2i, such that this relationship is expected to hold only as a fixed point 
relationship for ^ — > oo. Following the same argument as for generating the ansatz for 
the larger system, we may transform A''^^""' A'^' to A^'^^B''^^'^""^ Then the left-hand side of 
the fixed point relationship simplifies, using right normalization, to A^'^^A^'^^' = p^^, and it 
becomes p^^' - pj^ If this relationship holds to high accuracy, the thermodynamic fixed 
point has been reached. One way of measuring the closeness of the two density operators 
is given by the fidelity ll65ll 



F(p'^,pr')-Tr 



PlPa 



(342) 



Inserting the definitions and using cyclicity properties of the trace, one can show that 
F(p^[\p^^^^^) - 2/ where s, are the singular values of A^'^'a'^"'^"''. 

Of course the algorithm can be stopped at any time, but then we have a finite system result which 
definitely can be improved by using finite-system DMRG. The convergence criterion given really 
gives us access to the thermodynamic limit state, which we might write down formally as 

(343) 

where we take f to be the iteration step when the convergence criterion is met. The question 
is now how to evaluate expectation values. Obviously, we cannot write down a finite network 
contraction as before; it will be of infinite size and therefore cannot be contracted naively. A 
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contraction can only work if we can reduce the number to a finite number of contractions. For 
finite-system networks we saw that left- and right-orthonormality allow to eliminate most con- 
tractions: For observables on sites / and j, one only had to consider the contractions on and 
between these two sites. There is however no reason why the thermodynamic limit state should 
meet normalization criteria; in fact, usually it does not. We therefore need an orthonormalization 
procedure. After that, expectation values can be evaluated as for a finite lattice with left- and 
right-orthonormalization. Because this procedure is also important for the next algorithm and 
conceptually a bit more advanced, I postpone it to an extra section. 

At this point it should be mentioned that iDMRG can be related to earlier algor ithmic ap- 
proaches under the name of PWFRG (product wave function renormalization group) 1 138 , 139l 
1 140.1 which alr eady contain part of the above ideas; iDMRG takes them to their natural comple- 



tion 1 141, 1421 



10.4. iTEBD: time evolution in the thermodynamic limit 

In this section, I switch to the FA notation, although the formulae can be easily translated into 
the A, Z?-formulation. Using our state fragment At'^^'lr^'"''^ At^^Ff'"'? , we can set up an infinite 
chain 

(A'*F'^'^' A*r*'^'+' ){^js^^Y^<^i*2 j^B^Ba,,,-^ _ ^ _ 1^,^^ (344) 

(T 

just like in the previous section, where F'*'^ = T''^'^ , F*'^ = F^^^?, A"* = K^^-^^ and A* = K^^. 
The fragment may be the result of a converged ground state calculation or from some simple 
starting state that one can construct exactly. 

We can now write down a time evolution in the Trotter form by applying an infinitesimal time 
step to all odd bonds (which I will refer to as AB) and then on all even bonds (which I will refer 
to as BA). The bond evolution operators will be exactly as in the tMPS/tDMRG/TEBD cases, I 
will refer to them as Uab — Zlo-^o-so-^o-' 

U'2'p-'^'''^'>\o-AcrB){a-'^o-'g\ and similarly Uba- 
As we have already seen [cf. Eq. ( I338l l1. a full two-site fragment consists of a product of five 
matrices, A'^F'*""' A*F*°"* A'*. Then time evolution on bond AB yields a set of matrices 

l^<TAas ^ f/J^'^''''">«A^F'*'^^A*F^'"«A'*. (345) 

Upon the by now standard reshaping we obtain by S VD 

M'^'"^« = A'^''A*B'^«, (346) 

where the new A* (and correspondingly A"^' and B"'") are truncated as in tDMRG or TEBD, to 
replace the old one. Using A'* (still from the last iteration), we can define new F'^""' and F*""" 
(via A'^' = A'^F'*'^" and B'^* = F*'^* A^*). This defines a new "unit cell". 

If we write it down and attach another one, we can read off the bond BA in the center of the 
two AB unit cells as a^F^°"''A'^F'*°"'* A*, for which time evolution gives 

A?-«-^ = ^ f/^f -^'"^'^A^F^'^^A^F^^^A^. (347) 



Reshaping and SVD gives us 
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(348) 



where again A"* (and correspondingly the other matrices) are truncated and replace the old ones. 
Using A* (still from the last iteration), we can define new T*""" and T'^'^' (via B""" = T'*""" A* and 
A""" = A^r*"""). The problematic division by^small singular values can be avoided by the simple 
modification already discussed for TEBD 1 105]. 

By applying sequences of infinitesimal bond evolution operators we can therefore set up a 
real or imaginary time evolution for the thermodynamic limit. This algorithm is referred to as 
iTEBD, because it provides the infinite-size generalization of TEBD |63]. 

Again, the question of orthonormality arises 16411 . Let us assume that the initial state was 
meeting orthonormality criteria. A pure real-time evolution generates a sequence of unitaries 
acting on the state, which preserves orthonormality properties. But the inevitable truncation 
after each time step spoils this property, even though truncation may only be small. To turn 
this method into a viable algorithm, we have to address the issue of orthogonalization in the 
thermodynamic limit, as for iDMRG, after each step. 

Let me mention here that McCulloch t65il has shown that iDMRG can be turned into iTEBD 
by replacing the minimization on the central bond by a time-evolution on the central bond, with 
some conceptual advantages over the original iTEBD algorithm. 

Let me conclude this section by a few words on extrapolation. In finite-system DMRG (or 
MPS), the recipe was to extrapolate for each finite length L results in D to maximize accuracy, 
and then to extrapolate these results in L to the thermodynamic limit. Here, we are working di- 
rectly in the thermodynamic limit (assuming that iDMRG has been taken to convergence), and the 
extrapolation in D remains. Interestingly, at criticality, this extrapolation i n D h as profound and 
useful connections to entanglement entropy scaling and critical exponents Effectively, the 

finite matrix dimension introduces a finite correlation length into the critical system, not unlike 
the finite-system DMRG case, where the length scale on which the MPS -typica l supe rposition of 
exponentials mimicks a power-law properly also scales with some power of D i 14411 . 

70.5. Orthogonalization of thermodynamic limit states 

Within iDMRG on a finite system. A- and B-matrices retain left- and right-normalization; this 
implies that the left and right block states are orthogonal among each other, as shown previously. 
We will call a state with this property orthogonal in a slight abuse of conventional usage. As we 
have seen in the previous section, we may use a fragment 

^[fKA[nBm<r?(A[f-i])-i (349) 
that we can repeat to build up an infinitely long chain, 

= 2 . . . A'"'A[^1b'^'+'(A[^-i])-1A'"'-=AMb'^'-'(A"'"")"'A'"''-''AMb'"''-=(A[^-i])-i . . . |o-), 
cr 

(350) 

where I have simplified the notation of A, B. The problem with these states is that, for an 
arbitrary bipartition into two blocks, the states on the left and right blocks will in general not be 
orthogonal: if we transform A'^^^B into AJS^^ as described above, the chain will read 

= ^ . . . A'^'A^^f A'^-'A'^-'f A'^'A'^op . . . |cr), (351) 
cr 

where P = Kf{AS''-'^^)-\ If we absorb P into the A to its left, AP — > A, the normalization 
condition becomes 

^ T^'T = P"'' ^ A^^'A'^P = P' P. (352) 

(T cr 

HI 



In general, however, P^P + I. This is not only the case if ( is small and we are far from the 
infinite-system fixed point. It is also the case at the fixed point as long as the discarded state 
weight is finite, which is usually the casejn DMRG calculations, even if it is very small. 

As pointed out by Orus and Vidal|64] - in the presentation I follow fHs"] -, a condition to 
detect orthonormality is to check whether the expectation value of the unit operator between two 
block states \a), \a') is 6a.a' (see Fig.l62Ti. Let us consider an expectation value contraction as for 
a finite system and assume we have carried it out up to site 0, coming from the left, -oo. The 
result will be a matrix-like object Ca'a, corresponding to the open legs. Let us now consider the 
operation EiiC), which carries the contraction two steps further, i.e. over sites 1 and 2. This 
transfer operator reads [cf. Sec. 14.2.211 

£^(C) = ^ P^A''-^A''''^CA'^'A''-P. (353) 

0"|0"2 

For an orthonormal state, we want that El{I) = /, which is just the expectation value matrix 
the unit operator produces for orthonormal block states. What we get, however, is, using the 
left-normalization condition, El{I) = P' P. As the system extends to infinity, Ei(I) - I must be 
associated with the largest eigenvalue; normalizability of the entire state implies that the largest 
eigenvalue must be 1 . The "quadratic" form of E^ implies that the associated eigenmatrix Vl is 
hermitian and non-negative. An eigenvalue or singular value decomposition allows to decompose 
Vl - X^X, where X is invertible. We can insert X^^X after each P, such that the unit cell becomes 
XA°"' A°"2f X ' and the new transfer operator reads 

E'^iC) = X'^'P'A°''^A'^'^X^CXA'^'A'^'-PX-K (354) 

0-10-2 

Then 

£[(/) = X-^'VlX-^ = / (355) 

from the eigenmatrix properties of Vl with respect to El- (If the largest eigenvalue of El happens 
not to be I, X must be suitably scaled.) 

Inserting the definition of P in XA"'' A"'- PX^\ undoing the transformation to A, and trans- 
forming A'^'Al^l ^ A^[^B°'\ the unit cell becomes XA^'^B^'^B'^HA^'-'^Y'^X-^ Shifting the 
starting point of the unit cell it becomes (A^'^-^Y^X-^XA^l^ B°'' B"^'- = QB°''B'^\ where Q = 
(A^(-^^)-^X'^XA^l^ - (A^''-^^)-^A^l\ independent of X. Calculating a contraction from the right 
leads to a transfer operator 

Er(C) = Yj QB'^'B'^-CB'^'-'B'^'^Q'. (356) 
0-10-7 

The same eigenvalue/eigenmatrix argument as before leads to the dominant eigenmatrix Vr - 
YY\ Y invertible, and a unit cell Y-^QB'^'B'^^Y. This in turn leads to E'^^iQ with E'j^{I) = /. 

If we insert the definition of Q into the unit cell, return from to A°"' and make Q explicit, 
the unit cell reads Y'\A^^^^^)~^X~^XA°'' A^^^B°'-Y, shifting its origin we obtain 

XA'^'A^^^B'^'-YY-\A^^-^^)-^X-\ (357) 

which can be brought back to the original form of the unit cell by setting A""' <— XA"'' , B""- <— 
B'^-Y and A''^"'^ <— XA^^^^^Y, but now with proper left- and right-normalization ensured. 
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Figure 62: If matrices are properly normalized, the thermodynamic limit ansatz generates both a left- and right- 
normalization condition. Note that the two transfer operators are defined on two differing two-site unit cells of the 
thermodynamic limit state. 



More precisely, the new unit cell leads to Ei{I) - I and £«(/) - I. But note that Ei and 
Er are constructed from slightly shifted unit cells, namely A°"' A''^^B°"2(A'^"'^)"' for Ei and 
(^[^-i])-i^o-i^[^l/jo-2 for Er, as shown in the pictorial representation. We can lump together 
the first and second unit cells into left- and right-normalized two-site matrices A°"'°"' and Z?°"'°"-. 
These can now be decomposed into left- and right-normalized matrices in the standard way, giv- 
ing Af''°^'A'^^°"- and B[i]°"ifiPlo"2 Note that these matrices are of course different from those we 
had originally. 

The thermodynamic limit state can now be written as 

. . A^'J^^'API'^^A^'^'^mP)'^^ . . . Io-) (358) 

cr 

or analogously using B'^'^^, for left- and right-canonical form. Of particular interest for expecta- 
tion values is the mixed-canonical form, with A-matrices on the left and B-matrices on the right. 
If we consider the two underlying unit cells, we see that at the boundary, both want to incorporate 
the same (A'^"'^)"' to generate A and Z?-matrices. This problem can be solved by inserting the 
identity / = A'^'^"''(A'^'^"'1)"' after the problematic (A^^"'')"'. Then we can immediately write 
down a mixed-canonical form as 

|(/,) = ^ . . . A['l'^'Al2J'^=At'l'^MPJ'"-'A[''-'JB['J'^5BPKf,^[i]cr,^[2Ks ^ ^ ^ 1^^^ (359) 
cr 

10.5.1. Calculation of expectation values in the thermodynamic limit 

In finite systems, we have seen how an expectation value can be calculated by transferring 
an object C''', starting as a dummy scalar C'^"' = 1 from C^'^ to C^'^'^ by means of a transfer 
operator e'^^, where hatO is the locally acting operator in the expectation value structure (mostly 
the identity, except, say, two sites if we are looking at a two-point correlator). In the case of the 
identity operator, for left-normalized matrices, the transfer operator mapping from left to right 
maps the identity to the identity; similarly, the transfer operator mapping from right to left maps 
the identity to the identity if formed from right-normalized matrices. 

The same structure has been found in the last section for the thermodynamic limit state and 
its two two-site transfer operators Ei and Er. This allows a dkect transfer of the old results. 
Assume we want to calculate {il/\OiOi\^)\ then we bring the state into the mixed-canonical form 
of Eq. ( 13591 ), with A-matrices up to site / or / -i- 1 (depending on the odd-even structure), then 
A^^"'^, followed by B-matrices up to infinity. Contracting from the left over all A-matrices up 
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Figure 63: Exploiting left- and right-normalization, the evaluation of an infinite contraction can be reduced to a finite 
contraction that can be dealt with as for any finite MPS. 

to and from the right all B-matrices, we obtain a remaining finite network as in Fig.|63] This 
expectation value is then evaluated as in a finite network; if we use the C''' and transfer operator 
notation, it starts from C'"^ - I. The difference is that at the end, A'^ shows up in the contrac- 
tion and the final reduction to a scalar is done by the closing Saa' line (which can be read as a 
trace): assuming the last site is /, then the final expectation value is given in the last step as 



where we have used the relationship between A-matrices and reduced density operators in canon- 
ical representations. 

This calculation can be reduced easily to the case of the overlap of two states, which is just the 
matrix element of the unit operator between them. Assume that both states are in translationally 
invariant form, e.g. by using left-normalized matrices A^"^^ and A^^^ (and A^^^, A^^ respectively). 
We now carry forward an infinite overlap calculation by two sites (say 1 and 2) towards the right 
using El. If the current overlap matrix is C, it is carried forward as 



If we decompose C in the eigenmatrices of Ei, in the thermodynamic limit only the largest 
eigenvalue contribution will survive. For an orthonormal state, for the overlap with itself, C - I 
and A - \ are the dominant eigenpair A smaller A in the overlap of two states can be interpreted 
as an overlap per site, while of course the two states are orthogonal with respect to each other in 
the thermodynamic limit (overlap lim^^oo A^ - 0). Such thermodynamic overlaps (or fidelities) 
per site can be used very nicely to detect quantum phase transitions by overlapping ground states 
for two Hamiltonians with slightly different parameters L65.1 . 

11. Conclusion: other topics to watch 

After moving through this long list of topics, focussing on the fundamental algorithmic build- 
ing blocks, ground state searches, thermodynamic limit algorithms and a wealth of real and 
imaginary time methods at zero and finite temperature, the possibilities of DMRG and MPS- 
based algorithms are far from being exhausted. A few topics that I have not touched upon, 
but which I would like to mention briefly (again in a non-exhaustive list), are: transfer matrix 
DMRG methods (cf. the introductory section for references), DMRG and MPS with periodic 
boundary conditions, and as the most recent addition, MPS for continuous space ll66ll . which 
emerge as a beautiful generalization of coherent states and should allow for interesting applica- 
tions in field theories. For periodic boundary conditions quite a lot of results already exist, so 



{ilf\0' ® . . .OV) = TrAl^-'J-' cI'lAl^-'l = TrA^'-'^A^'-'^^C^'^ = Trp^C^'l 



(360) 




(361) 
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let me give just a brief overview. PBC have already been treated in the DMRG framework by 
intro ducing one long-ran ged interaction between sites 1 and L on an open-boundary chain (see 



e.g. 111451 Il46lll47l 114811 : however, the scaling of accuracy was consistently found to be much 
worse than for open boundary conditions. The underlying reason is (roughly speaking) that on 
a ring the surface between A and B doubles, hence the entanglement; given the exponential re- 
lationship to the MPS dimension, this means that resources have to go up from D to up to D^, 
meaning that for similar accuracy, the algorithm needs the square of time (sometimes referred 
to as D^-scaling, referring to the open boundary condition D). The physically adequate ansatz 
for MPS for periodic boundary conditions is given by Eq. (fTTT l: one needs roughly the same D 
as for OBC, but rerunning the variational ground state search algorithm on it scales as (be- 
cause the simplification of vectors instead of matrices on sites 1 and L does not occur) 1.56,1 . At 
the same time, the simplification of the generalized to a standard eigenvalue problem does not 
occur, which may lead to bad conditioning. A nice feature of the MPS representation for PBC 
is that one can generate eigenstates of momentum: For k - n{2n/L) and a (non-translationally 
invariant) MPS lij/} = Yjcr Tr (A^'|° "' . . .A^^^°'^)\(t}, the following state is a translationally invariant 
eigenstate of momentum k: 114911 

L-l 

^ e'*"Tr(A['J'^'+" . . . At^l'^^+'Ok)- (362) 



11=0 



Recently, interesting proposals to improve the scahng have been made |j57, 58], and this is a 
field of ongoing interest. Reference ll69ll discusses this topic quite extensively. 

1 think one may conclude by saying that while the fundamental framework of MPS is by now 
very well established, and while DMRG has come of age as one of the most powerful numerical 
methods available for strongly correlated quantum systems, even in the well-established field of 
one-dimensional systems many of the algorithms presented will still allow further improvement, 
bringing new applications into our reach. It is in fact quite surprising that for quite a few of the 
methods presented (and also the others) very little is known about their detailed behaviour in real- 
world problems, analyzing which might give interesting further ideas. Also, the ratio between 
applications done and applications doable seems very favourable for future exciting research. 
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